clear all
close all

% =====================================================================
% ====  "Risk-Sharing in Village Economies Revisited", by Bold and Broer

% ====   This programme solves the limited commitment village economy for
% both the standard specification (individual deviations) and the
% alternative specification (coalitional deviations)

% ====    Written by Tessa Bold and Tobias Broer

% ====    This version June 2021
% =====================================================================


filetosave='Vilcoal_24_10_20_vill_cd_hetinc'; % name of workspace saved

%====================================================
% select parameters
% ===================================================

village=1;                         % village indicator: 1 - Aurepalle; 2-Kanzara; 3-Shirapur
coaldev=1;                         % set to 1 for coalitional deviations, to 2 for self-insurance (SI), to 0 for individual deviations (ID), to 3 for full-insurance with coalitional deviations
Tessa1Tobi2=1;                     % sets paths depending on user below
unconditional=0;                   % set this to 1 for unconditional income data, to 0 for income data conditional on time fixed effets / demeaned data period-by-period
fixedeffect=1;                     % set this to 1 for an income process with fixed effects
heterogeneity=1;                   % set to 1 if you want to allow for heterogeneity in income processes of villagers
one_per_aut=0;                     % set this to 1 if you want deviations to imply one period of autarky
ytaxrate=0;                        % linear tax on income
scatter = 0;		           % set to 1 to just draw a scatter plot
Rov_sb = 0;                          % set this to 1 if you want the outside option for the Rov in the case of coalitional deviations to equal the average value from the solution to the next smaller village (second-best), rather than
                                   % just the utility from consuming average ROV income (first-best)
estrho1=0;			   % set this to 1 if you would like to estimate persistence, default is 0
rhogrid=0;                         % set to 1 if you want to have also a grid for CRRA rho - only to show non-identification in Aurepalle (village=1)
Nphi=1;                            % number of borrowing limits in SI economy - not used currently
server=0;                          % set to 1 to run on server; o wise path for workstations is used
buildup=0;                         % set to 1  to calculate ID model for subset of group sizes, not just whole village
Nincproc=1;                       % number of income processes to look at
if estrho1==0
maxincerr=2/3;                     % maximum income measurment error in multiples of income variance (in levels)
elseif estrho1==1
maxincerr=0;
end
model_var_dy=1;                    % set to 1 if you want to use the model-implied variance of dlog(y) to specify maximum meas error, otherwise var(dlog(y)) from data is used, slightly different as we target var(log(y)), not var(dlog(y())
T=6200; Tinit=201;                 % number of simulation periods in total, T-Tinit+1 equals number of initial periods to discard
maxcerr=0.9;                       % max cons meas error in log levels
if scatter==0
N_cerr=0;                         % number of support points for cons meas error
elseif scatter==1
N_cerr=0;
end
if coaldev==0
    convlength=0;                  % the conv criterion must be fulfilled at least in convlength consecutive periods to avoid jumping 'in'
else
    convlength=15;
end
if village==1
    vss_high_inddev=34-1;          % the number of total village members is vss_high_inddev+1 under individual deviations
    vss_high_coaldev=34-1;
elseif village==2;
    vss_high_inddev=37-1;                 
    vss_high_coaldev=37-1;
elseif village==3
    vss_high_inddev=31-1;               
    vss_high_coaldev=31-1;
end
vssbuildupvec_inddev=[1:7,10,15,vss_high_inddev]; % for buildup==1, this is the vector of village sizes for which the solution is calculated in the individual deviations case
if server==1

else
    if Tessa1Tobi2==2
        outputpath=sprintf('C:/Users/tbroer/Dropbox/Research/Current_Projects/tessa/JEEA Replication/Programmes/Output estimation August 2018 conditional')
        programpath=sprintf('C:/Users/tbroer/Dropbox/Research/Current_Projects/tessa/JEEA Replication/Programmes');
    elseif Tessa1Tobi2==1
        outputpath=sprintf('C:/Users/tbold/Dropbox/Tessa and Tobi/JEEA Replication/Programmes/Main Results and Robustness')
        programpath=sprintf('C:/Users/tbold/Dropbox/Tessa and Tobi/JEEA Replication/Programmes')
    end
end
        cd(programpath)
if buildup==0
    if rhogrid==0 && scatter == 0
        betavec0=[0.5:0.04:0.78,0.8:0.01:0.98];
elseif scatter ==1
betavec0=[0.83];
    elseif rhogrid==1
        betavec0=[0.83,0.85,0.87,0.89,0.91,0.93,0.95,0.97];
    end
elseif buildup==1

betavec0=[0.86,0.905,0.91,0.94]; %betavec0=[0.7:0.04:0.82,0.84:0.02:0.98];

end
if rhogrid==0
if scatter==1
betavec1=[0.94];
else
 betavec1=[0.84:0.02:0.88,0.89:0.01:0.96,0.9625:0.0025:0.99];
 betavec1=[0.9:0.01:0.98];   %%this is your betavec for income heterogeneity
end
elseif rhogrid==1
    betavec1=[0.84,0.86,0.88,0.90,0.92,0.94,0.96,0.98];
end
betavec3=[0.96:0.0025:0.99];
if coaldev<2 || coaldev==3
    
    rhovec=[1];
elseif coaldev==2
    rhovec=ones(1,Nphi);
end

if rhogrid==0
    gapcritset=0.0001;              %convergence criterion for change in values
elseif rhogrid==1
    gapcritset=0.001;
end

rov_average=1;                      % set to 1 to specify the rest of village in terms of average
                                    % consumption
momentclevinlog=1;                  % set this to 1 if you want to calculate the moments for levels of
                                    % consumption based on log consumption
maxitgapset=400;                    % maximum number of iterations
Pun_fac=0.0;                        % consumption penalty in the first period after deviation equal
                                    % to caut*Pun_fac
% grid of time-varying planner weights
lambda=[0.05:0.01:0.95]';
m=length(lambda);
if floor((m+1)/2)~=((m+1)/2)
    disp('care: lambda has to have uneven no of points in first dimension')
    stop
end
%S_KEEP=nan(100*3,2,vss_high_inddev,4,3);               % pre-allocate the matrices for coalitional deviations, so that you don't run simulations vss times
S_rov_KEEP=nan(100,1,vss_high_inddev,4,3);              % this one is for a maximum of vss=19 for the CD
P_rov_KEEP=nan(100,100,vss_high_inddev,4,3);

% ====================================================
% Specify Income process for individual and Rest of Village (rov)
% ====================================================
% data moments for log income from ICRISAT data - all villages
% first col: Aurepalle, second: Kanzara, third: Shirapur
if fixedeffect==1
    if unconditional==0 & heterogeneity==0
    datafile='armoments_fe0';
    elseif unconditional==0 & heterogeneity==1
    datafile='armoments_het_fe0';   
    elseif unconditional==1
    datafile='armoments_fe1';
    end
elseif fixedeffect==0
     datafile='armoments';
end
if server==0 & Tessa1Tobi2==2
    cd('C:\Users\tbroer\Dropbox\Research\Current_Projects\tessa\JEEA Replication\Data\Output')
elseif server==0 & Tessa1Tobi2==1
    cd('C:\Users\tbold\Dropbox\Tessa and Tobi\JEEA Replication\Data\Output')
end
eval(['load ' datafile '.out;']);
eval(['covyylag=' datafile '(1,:);']);
eval(['vardy=' datafile '(2,:);']);
eval(['vary=' datafile '(3,:);']);
eval(['meany=' datafile '(4,:);']);

if model_var_dy==1
    vardy=[0.231,0.125, 0.269];
end
cd(programpath)
if maxincerr>0
if Nincproc>1
varnu_max=maxincerr*vardy/2;                % nu is measurement error in log levels
else
  varnu_max=zeros(1,2);    %%%you're doing heterogeneity in incomes now, for one village   
end

rho_min=covyylag./vary;                      % AR(1) coefficient w/o meas error
rho_max=covyylag./(vary-varnu_max);          % AR(1) coefficient w max meas error
for g=1:2
rho1vec(g)=linspace(rho_min(g),rho_max(g),Nincproc);
end
varnu=max(vary-covyylag./rho1vec,0); % max prevents the first entry to be -e17
vareps=(vary-varnu).*(1-rho1vec.^2);
else
    rho1_data=covyylag(village)/vary(village);
    vareps_data=vary(village).*(1-rho1_data.^2);
    rho1vec=[-0.6,-0.4,-0.2,0,0.2,0.4,0.6];
    vareps=vary(village).*(1-rho1vec.^2);
    varnu=zeros(1,length(vareps));                % nu is measurement error in log levels
    Nincproc=length(vareps);  
end

%%%we'll have gmatch and kmatch, for kmatch, the loop will have to simulate
%%%for both kmatch, but separate simulations for each gmatch
matching=[1 1 2 2; 1 2 2 1; 3 3 4 4; 3 4 4 3];

% ===================================
% Iterate over income processes
% ===================================
    disp('now solving for village,coaldev, unconditional, Nincproc, fixedeffect one_per_aut Rov_sb - press any key to continue')
    disp([village coaldev unconditional Nincproc fixedeffect one_per_aut Rov_sb])
    pause
indicator1=zeros(1,size(rhovec,2));
incproc=1;
    if rhogrid==0 && buildup==0 && maxincerr>0 && Rov_sb==0 && one_per_aut==1 && scatter == 0 
        filetosaveest=['Momentsest' int2str(village) int2str(coaldev) int2str(incproc) int2str(unconditional) int2str(fixedeffect) int2str(one_per_aut)];
    elseif maxincerr==0
        filetosaveest=['EstrhoMomentsest' int2str(village) int2str(coaldev) int2str(incproc) int2str(unconditional) int2str(fixedeffect) int2str(one_per_aut)];
elseif one_per_aut==0 % these are the benchmark options
        filetosaveest=['Inc_het_Momentsest'];
    elseif Rov_sb==1
        filetosaveest=['SBMomentsest' int2str(village) int2str(coaldev) int2str(incproc) int2str(unconditional) int2str(fixedeffect) int2str(one_per_aut)];
    elseif rhogrid==1
        filetosaveest=['GridMomentsest' int2str(village) int2str(coaldev) int2str(incproc) int2str(unconditional) int2str(fixedeffect) int2str(one_per_aut) ];
    elseif buildup==1
        filetosaveest=['BuildupMomentsest' int2str(village) int2str(coaldev) int2str(incproc) int2str(unconditional) int2str(fixedeffect) int2str(one_per_aut)];
elseif scatter==1
 filetosaveest=['ScatterMomentsest' int2str(village) int2str(coaldev) int2str(incproc) int2str(unconditional) int2str(fixedeffect) int2str(one_per_aut)];

    end

    % set minimum and maximum village size
    if coaldev==1                       
        if rhogrid==0
            vss_set=[1:5];% prior solutions have shown that 1-7 is the relevant set. Otherwise use, for example, [1:12,17,22,27,vss_high_coaldev]
        elseif rhogrid==1
            vss_set=[1:5]
        end
        vss_high=vss_set(end);
    elseif coaldev==0
        if buildup==0
            vss_set=[vss_high_inddev];
            lagind=1;
            vss_high=vss_set(end);
        elseif buildup==1
            vss_set=vssbuildupvec_inddev;
            lagind=1;
            vss_high=vss_set(end);
        end
    elseif coaldev==2;
        vss_set=vss_high_inddev;
        vss_high=vss_set(end);
    elseif coaldev==3
        vss_set=[1:10];
        vss_high=10;
    end
    
    
    % ===================================
    % iterate over discount factors
    % ===================================
    if coaldev==0
        betavec=betavec0;
    elseif coaldev==1
        betavec=betavec1;
    elseif coaldev==2
        r=0.04;
if scatter==0
        betamax=1/(1+r)-0.015;
        betavec=linspace(0.8,betamax,20);
elseif scatter ==1
betavec=0.92;
end
        Z=incomevals;
        phinat=min(Z)/(r-1);
        if Nphi>1
            phivec=linspace(phinat,0,Nphi);
        else
            phivec=0;
        end
    
    end
    N_cerr=0;
    Momentsest=nan(35,N_cerr+1,size(betavec,2),1,vss_high_inddev);
    for Qbeta=1:size(betavec,2)
        beta=betavec(Qbeta);
        if rhogrid==1
            if coaldev==0
                if beta==0.83
                    rhovec=1.468461538461538
                elseif beta==0.85
                    rhovec=1.316666666666667
                elseif beta==0.87
                    rhovec=1.16
                elseif beta==0.89
                    rhovec=1
                elseif beta==0.91
                    rhovec=0.8340
                elseif beta==0.93
                    rhovec=0.6611
                elseif beta==0.95
                    rhovec=0.4808
                elseif beta==0.97
                    rhovec=0.2912
                end    
            elseif coaldev==1
                if beta == 0.98
                    rhovec=0.366222222222222
                elseif beta==0.84
                    rhovec=2.389
                elseif beta==0.86
                    rhovec=2.11
                elseif beta==0.88
                    rhovec=1.821
                elseif beta==0.90
                    rhovec=1.5515
                elseif beta==0.92
                    rhovec=1.28718
                elseif beta==0.94
                    rhovec=1
                elseif beta==0.96
                    rhovec=0.6935
                else
                end
            end
        end

        % ===================================
        % iterate over CRRA
        % ===================================
        for Qrho=1:size(rhovec,2)
            rho=rhovec(Qrho)  ;
            if coaldev==1
                lagind=[];
            else
                lagind=1;
            end
            
            
           
            STABLEHET=[];
            gmatch=1 
                if gmatch==1
                    gmatchindex=1
                end
                P_LAG=nan(1000,1000,length(vss_set));Y_lag=nan(1000,1000,length(vss_set));
                waut_lag=[]; waut_avg_lag=[]; waut_rov_sb=[]; waut_rov=[];waut_temp=[]; 
                
                   
                % specify a markov process for individual income
                for k=[0 1]
                [incomevals, Q1] = rouwenhorst(rho1vec(gmatchindex+k),vareps(gmatchindex+k)^0.5,3);   
                incomevals=exp(incomevals');
                incomevals_match(:,k+1)=incomevals;
                Q1_match(:,:,k+1)=Q1;
                end
                for k=1:2
                QQQ_match(:,:,k)=Q1_match(:,:,k)^1000;
                SS_match(:,k)=QQQ_match(1,:,k)';
                incomevals_match(:,k)=incomevals_match(:,k)/(SS_match(:,k)'*incomevals_match(:,k));
                end
                %%%adjustment: Multiply income process of the rich by rich/average income, of the poor by poor/average income. Make sure that overall it sums to 1 across the stationary distribution.
                for k=1:2
                incomevals_match(:,k)=incomevals_match(:,k)*meany(gmatchindex+(k-1))/mean(meany);
                end
                for k=1:2
                cumQ_match(:,:,k)=cumsum(Q1_match(:,:,k),2);
                end
                
                
                
                incomevals_match=(1-ytaxrate)*incomevals_match+ytaxrate; % tax as in Krueger and Perri JET 2011
              
                
                              
             STABLE=[];
             if Qbeta<=5
                 vss_set=[1 2 3 4];
             else
                 vss_set=[1 2 3 4 5 6];
             end
            for vsscount=1:length(vss_set)       % iterate over village size
                
                vss=vss_set(vsscount);
                vss_max=max(vss_set);
                if vsscount>1 & coaldev==1
                    lagind=lagind+(vss-(vss_set(vsscount-1)+1));
                end
                disp('now solving for parameters')
                disp([Qbeta,Qrho])
                disp('vss')
                disp(vss)
                
                
                    % this is the scaling factor: if you want RoV consumption as average consumption per HH, scale rov variables
                    % by vss
                    if rov_average==1
                        ScaleRov=vss;
                    else
                        ScaleRov=1;
                    end
                    
                               
                    MATCHES=nmultichoosek([1 2],vss);                    %%%who does the individual match with in the rest of the village
                    MATCHES_LENGTH(vss)=size(MATCHES,1);
                    MATCHESKEEP(1:vss+1,1:vss,vss)=MATCHES;
                    for kmatch=1:2  %%%whether the individual has the `low' or the `high' value 
                    for villmatch=1:length(MATCHES)
                    
                    n_type(1)=length(find(MATCHES(villmatch,:)==1)); %%%how many of each income type in rov
                    n_type(2)=length(find(MATCHES(villmatch,:)==2)); %%%%how many of each income type in rov    
                       
                                          
                    csim=[];  cc =[]; eeewaut =[]; eewaut =[]; ewaut =[]; waut =[]; waut_check=[]; waut_check_ind=[]; rr =[]; cshare=[]; raw =[]; cshare =[]; dcshare =[]; dcshare_raw =[]; S_rov=[]; P_rov =[];
                    S_rov_new =[]; P_rovnew =[]; S =[]; P =[]; caut =[]; c =[]; waut =[]; gammagrid =[]; gammar =[]; phi =[]; phisol =[]; csol =[];csolold=[]; cfirstbest =[];
                    wfirstbest =[]; ewfirstbest =[]; wsbnewsol=[];
                    Y =[]; vs=[];  idio_dist_autnew=[];idio_dist_aut_new=[];normperf=[];w1=[];indSj=[]; indSjj=[]; Ssim=[]; aggvillincome=[];aggincgrowth=[];aggvillincgrowth=[];
                    S_rov_aut=[]; S_rov_sb=[]; S_rov_aut_all=[];S_rov_aut_all2=[];
                    w =[];  phinew =[]; Yagg =[]; wint =[]; wsb =[]; ewsb =[]; csol =[]; phisol =[]; nonconv =[];   gammarfun =[];  x0 =[]; index =[]; ewsb =[]; bind=[];
                    waut=[];
                    
                    % ===================================
                    % build the state space of ROV income
                    % ===================================
                    idio_dist=[];
                    idio_dist_aut=[];
                    idio_dist_aut_all=[];idio_dist_aut_all2=[];
                    idio_dist_sb=[];
              %%%For income heterogeneity it is like the rest of the village consists of two types of households we need to keep track of       
                       
                       if n_type(1)>0
                       for k1=0:n_type(1)
                        for k2=0:n_type(1)
                            for k3=0:n_type(1)
                                II=zeros(1,3);
                                if k1+k2+k3==n_type(1)
                                    % the vector of aggregate incomes from
                                    % permutations of vss households
                                    S_rov=[S_rov; k1*incomevals_match(1,1) + k2*incomevals_match(2,1) + k3*incomevals_match(3,1)];
                                    % the corresponding distribution of
                                    % individual income values
                                    idio_dist=[idio_dist; k1 k2 k3];
                                    
                                    
                                    
                                    end
                                end
                            end
                       end
                       S_rov_1=S_rov;
                       idio_dist_1=idio_dist;
                       end
                       S_rov=[];
                       idio_dist=[];
                       if n_type(2)>0
                       for k1=0:n_type(2)
                        for k2=0:n_type(2)
                            for k3=0:n_type(2)
                                II=zeros(1,3);
                                if k1+k2+k3==n_type(2)
                                    % the vector of aggregate incomes from
                                    % permutations of vss households
                                    S_rov=[S_rov; k1*incomevals_match(1,2) + k2*incomevals_match(2,2) + k3*incomevals_match(3,2)];
                                    % the corresponding distribution of
                                    % individual income values
                                    idio_dist=[idio_dist; k1 k2 k3];
                                    
                                    end
                                end
                            end
                       end
                       S_rov_2=S_rov;
                       idio_dist_2=idio_dist;
                       end
                       
                         
                       
                    
                    %%%combine S_rov_1 and S_rov_2
                    if n_type(1)>0 & n_type(2)>0
                    S_rov=[kron(S_rov_1,ones(size(S_rov_2))), kron(ones(size(S_rov_1)),S_rov_2)];
                    idio_dist=[kron(idio_dist_1,ones(length(idio_dist_2),1)), kron(ones(length(idio_dist_1),1),idio_dist_2)];
                    elseif n_type(1)>0 & n_type(2)==0
                    S_rov=S_rov_1;
                    idio_dist=idio_dist_1;
                    elseif n_type(2)>0 & n_type(1)==0
                    S_rov=S_rov_2;
                    idio_dist=idio_dist_2;
                    end
                     
                       
                    % resort S_rov in ascending order
                    [S_rov,S_rovind]=sortrows(S_rov);
                    idio_dist=idio_dist(S_rovind,:);
                    S_rov=sum(S_rov,2);
                    
                    
                   if vss>1
                        if n_type(1)>0 & n_type(2)>0
                        for k=1:vss
                            for g=1:length(idio_dist)
                                    index=0;
                                    for j1=0:vss-k
                                        for j2=0:vss-k
                                            for j3=0:vss-k
                                                for j4=0:vss-k
                                                    for j5=0:vss-k
                                                        for j6=0:vss-k
                                                if j1+j2+j3+j4+j5+j6==vss-k && (idio_dist(g,1)>=j1 & idio_dist(g,2)>=j2 & idio_dist(g,3)>=j3 & idio_dist(g,4)>=j4 & idio_dist(g,5)>=j5 & idio_dist(g,6)>=j6)
                                                    index=index+1;
                                                    temp=[j1 j2 j3 j4 j5 j6];
                                                    idio_dist_aut_all(g,:,index,k)=[idio_dist(g,:)-temp];
                                                    S_rov_aut_all(g,1,index,k)=[(idio_dist(g,1:3)-temp(1:3))*incomevals_match(:,1)]+[(idio_dist(g,4:6)-temp(4:6))*incomevals_match(:,2)];
                                                end
                                            end
                                        end
                                    end
                                end
                            end
                            
                        end
                                
                            end
                        end
                        elseif n_type(1)==0 || n_type(2)==0
                         for k=1:vss
                            for g=1:length(idio_dist)
                                    index=0;
                                    for j1=0:vss-k
                                        for j2=0:vss-k
                                            for j3=0:vss-k
                                                if j1+j2+j3==vss-k && (idio_dist(g,1)>=j1 & idio_dist(g,2)>=j2 & idio_dist(g,3)>=j3)
                                                    index=index+1;
                                                    temp=[j1 j2 j3 ];
                                                    idio_dist_aut_all(g,:,index,k)=[idio_dist(g,:)-temp];
                                                    if n_type(1)>0
                                                    S_rov_aut_all(g,1,index,k)=[(idio_dist(g,1:3)-temp(1:3))*incomevals_match(:,1)];
                                                    elseif n_type(2)>0
                                                    S_rov_aut_all(g,1,index,k)=[(idio_dist(g,1:3)-temp(1:3))*incomevals_match(:,2)];    
                                                    end
                                                end
                                            end
                                        end
                                     end
                                end
                            end
                        end
                    end
                    
                    
                    % ===================================
                    % simulate aggregate income transitions
                    % ===================================
                    %%%
                    tic
                    T_trans=50000;
                    aggincind=zeros(length(S_rov),length(S_rov));
                    for jjj=1:length(S_rov)
                        jjj;
                        if n_type(1)>0 & n_type(2)>0
                        inccount=zeros(T_trans,vss,6);
                        agginc=zeros(T_trans,6);
                        else
                        inccount=zeros(T_trans,vss,3);
                        agginc=zeros(T_trans,3);
                        end
                        for iii=find(idio_dist(jjj,:)>0)
                            if iii<=3 
                            if n_type(1)>0%%%this is the first `type' of individual in the rov
                            for kkk=1:idio_dist(jjj,iii)
                                rng(jjj*100+iii*10+kkk,'twister');
                                shock=rand(T_trans,1);
                                rr=(find(shock<=cumQ_match(iii,1,1)));
                                inccount(rr,kkk,iii)=1;
                                rr=(find(cumQ_match(iii,1,1)<shock & cumQ_match(iii,2,1)>=shock));
                                inccount(rr,kkk,iii)=2;
                                rr=(find(cumQ_match(iii,2,1)<shock & cumQ_match(iii,3,1)>=shock));
                                inccount(rr,kkk,iii)=3;
                            end
                            incomevals_aux=incomevals_match(:,1);
                                elseif n_type(2)>0 & n_type(1)==0 %this is when we only have type two individuals in ROV
                            for kkk=1:idio_dist(jjj,iii)
                                rng(jjj*100+iii*10+kkk,'twister');
                                shock=rand(T_trans,1);
                                rr=(find(shock<=cumQ_match(iii,1,2)));
                                inccount(rr,kkk,iii)=1;
                                rr=(find(cumQ_match(iii,1,2)<shock & cumQ_match(iii,2,2)>=shock));
                                inccount(rr,kkk,iii)=2;
                                rr=(find(cumQ_match(iii,2,2)<shock & cumQ_match(iii,3,2)>=shock));
                                inccount(rr,kkk,iii)=3;
                            end
                            incomevals_aux=incomevals_match(:,2);
                                end
                            agginc(:,iii)=sum(incomevals_aux(inccount(:,1:idio_dist(jjj,iii),iii)),2);
                            clear incomevals_aux;
                            elseif iii>3 %%%this is the second `type' of individual in the rov
                            for kkk=1:idio_dist(jjj,iii)
                                rng(jjj*100+iii*10+kkk,'twister');
                                shock=rand(T_trans,1);
                                rr=(find(shock<=cumQ_match(iii-3,1,2)));
                                inccount(rr,kkk,iii)=1;
                                rr=(find(cumQ_match(iii-3,1,2)<shock & cumQ_match(iii-3,2,2)>=shock));
                                inccount(rr,kkk,iii)=2;
                                rr=(find(cumQ_match(iii-3,2,2)<shock & cumQ_match(iii-3,3,2)>=shock));
                                inccount(rr,kkk,iii)=3;
                            end
                            incomevals_aux=incomevals_match(:,2);
                            agginc(:,iii)=sum(incomevals_aux(inccount(:,1:idio_dist(jjj,iii),iii)),2);
                            clear incomevals_aux;
                            end
                        end
                        agginc=sum(agginc,2);
                        for jjj2=1:length(S_rov)
                            qqq=(agginc<=S_rov(jjj2)+0.0000000001).*(agginc>=S_rov(jjj2)-0.0000000001);
                            aggincind(jjj,jjj2)=sum(qqq);
                        end
                        P_rov(jjj,:)=aggincind(jjj,:)/sum(aggincind(jjj,:));
                    end
                    disp('Simulation for agg income process takes')
                    toc

                    
                    if vss>1 && coaldev<2
                        for g=vss_set(find(vss-1>=vss_set))
                            
                                S_rov_aut_all(:,:,:,g)=S_rov_aut_all(:,:,:,g)/(g);
                            
                        end
                    end
                    S_rov=S_rov/ScaleRov;
                    
                    % ===================================
                    % state space: combinations of indiv and village income
                    % ===================================
                    
                    S=[kron(incomevals_match(:,kmatch),ones(size(S_rov))),kron(ones(size(incomevals_match(:,kmatch))),S_rov)];
                    idio_dist_rov=idio_dist;
                    
                    if vss>1 && coaldev<2
                        for g=vss_set(find(vss-1>=vss_set))
                             
                                for k=1:size(S_rov_aut_all,3)
                                    S_rov_aut_all2(:,:,k,g)=[kron(ones(size(incomevals_match(:,kmatch))),S_rov_aut_all(:,:,k,g))];
                                    idio_dist_aut_all2(:,:,k,g)=[kron(ones(size(incomevals_match(:,kmatch))),idio_dist_aut_all(:,:,k,g))];
                                   
                                end
                            
                        end
                        S_rov_aut_all=S_rov_aut_all2;
                        idio_dist_aut_all=idio_dist_aut_all2;
                        clear S_rov_aut_all2;
                        clear idio_dist_aut_all2;
                        clear S_idio_dist_aut_all2;
                    end
                    
                    P=kron(Q1_match(:,:,kmatch),P_rov);
                    
                    Y=S(:,1)+ScaleRov*S(:,2);			% possible aggregate income outcome
                    state=length(S);                    % number of states
                    agent=2;                            % number of 'agents' in the HH vs. RoV approximation equals 2
                    STATE(kmatch,villmatch,vss)=state;
                    
                    
                    % ===================================
                    % Compute set of autarky values
                    % ===================================
                    for k=1:agent
                        caut(:,k)=S(:,k);
                    end
                    waut=zeros(state,agent);
                    waut_avg_lag=zeros(state,agent);
                    waut_sb=zeros(state,agent);
                    % this is just the PDV of utility from consumption of
                    % income caut
                    eeewaut(:,1)= inv(diag(ones(1,state),0)-beta.*P)*utility(rho,caut(:,1));
                    eeewaut(:,2)= inv(diag(ones(1,state),0)-beta.*P)*utility(rho,caut(:,2));
                    
                    
                     gammagrid=[lambda vss/ScaleRov*(1-lambda)];
                        gammagrid=[gammagrid ones(m,1) gammagrid(:,2)./gammagrid(:,1)];
                        gammagridold=gammagrid;
                        clear gammagrid
                        gammagrid=[linspace(min(S(:,1)),1.2*max(S(:,1)),m)',linspace(1.2*max(S(:,2)),min(S(:,2)),m)',ones(m,1),linspace(1.2*max(S(:,2)),min(S(:,2)),m)'./linspace(min(S(:,1)),1.2*max(S(:,1)),m)'];
                        
                        GAMMAMIN(kmatch,villmatch,vss)=min(gammagrid(:,4));
                        GAMMAMAX(kmatch,villmatch,vss)=max(gammagrid(:,4));
                        
                        % calculate the autarky value for ROV from average utility
                        % of its members in village of size vss
                      

                        % in first iteration, or for individual deviations, waut is just consumption of income
                        if vss==1 || coaldev==0
                            for j=1:state
                                waut(j,1)=eeewaut(j,1);
                                % autarky values for the rest of village are always the same - perfect
                               	% risk-sharing among the rest of village
                                waut(j,2)=eeewaut(j,2);
                                	   

                            
                            end
                        else
                            % this calculates autarky values for the coalitional
                            % deviations case
                                                       
                            for j=1:state
                                waut(j,2)=eeewaut(j,2); %%%this is the actual outside option for the village
                                
                                %%%%%Here comes the big code that cycles through all
                                %%%%%potential coalitions of people, 
                                for g=1:size(S_rov_aut_all,4)
                                    
                                        for k=1:size(S_rov_aut_all,3)
                                            n_type_sub=[];
                                            if n_type(1)>0 & n_type(2)==0
                                            n_type_sub(1)=sum(idio_dist_aut_all(j,1:3,k,g));
                                            n_type_sub(2)=0;
                                            end
                                            if n_type(1)>0 & n_type(2)>0
                                            n_type_sub(1)=sum(idio_dist_aut_all(j,1:3,k,g));
                                            n_type_sub(2)=sum(idio_dist_aut_all(j,4:6,k,g));
                                            end
                                            if n_type(1)==0 & n_type(2)>0
                                            n_type_sub(1)=0;
                                            n_type_sub(2)=sum(idio_dist_aut_all(j,1:3,k,g));
                                            end
                                            gg=0;
                                            while gg<=g
                                                gg=gg+1;
                                                rr=find(MATCHESKEEP(gg,:,g)==1);
                                                if length(rr)>0
                                                n_type_matches(1)=length(rr);
                                                else
                                                n_type_matches(1)=0;
                                                end
                                                rr=find(MATCHESKEEP(gg,:,g)==2);
                                                if length(rr)>0
                                                n_type_matches(2)=length(rr);
                                                else
                                                n_type_matches(2)=0;
                                                end
                                                if n_type_sub==n_type_matches
                                                    villmatch_sub=gg;
                                                    gg=g+5;
                                                end
                                                clear n_type_matches
                                            end
                                            if STABLEHETFINAL(kmatch,villmatch_sub,g)==1
                                            S_lag_2=(S_rov_aut_all(j,1,k,g)==S_KEEP(:,2,kmatch,villmatch_sub,g)) ;
                                            S_lag_1=(S(j,1)==S_KEEP(:,1,kmatch,villmatch_sub,g));
                                            S_lag_ind=find(S_lag_1.*S_lag_2==1);
                                           
                                            agentstate=[];
                                            
                                            for kk=1:size(idio_dist_aut_all,2)
                                                if idio_dist_aut_all(j,kk,k,g)>0
                                                    agentstate=[agentstate kk];
                                                end
                                            end
                                            
                                            if size(S_lag_ind)>0
                                                for kk=1:length(agentstate) %%%note, you only need one entry per income state, not one entry per agent
                                                    
                                                        waut_check_ind(j,1,k,g)=WAUT_LAG(S_lag_ind,1,kmatch,villmatch_sub,g);
                                                        waut_check(j,kk,k,g)=WAUT_LAG(S_lag_ind,2,kmatch,villmatch_sub,g);
                                                    
                                                    
                                                end
                                            else
                                                waut_check_ind(j,1,k,g)=NaN;
                                            end
                                            else 
                                            waut_check_ind(j,1,k,g)=NaN;
                                            
                                            end
                                        
                                        end
                                    
                                end
                            waut(j,1)=max(max(max(waut_check_ind(j,1,:,:))));  
                            if isnan(waut(j,1))
                            waut(j,1)=eeewaut(j,1);
                            end
                            
                            end    
                        end
                        
                        % ====================================================
                        % Compute constrained insurance
                        % ====================================================
                        
                         
                       
                        % ====== make grids of weights, consumption, and Lagrange multipliers for updating ======
                        n=length(gammagrid);
                        w=zeros(n,state,agent);
                        gammar=zeros(n,state,agent);
                        c=zeros(n,state,agent);
                        phi=zeros(n,state,agent);
                        
                        % 'first-best' consumption given aggregate resources and planner weights
                        for j=1:state
                            c(:,j,1)=Y(j)./(1+ScaleRov*gammagrid(:,4));
                            c(:,j,2)=(Y(j)-c(:,j,1))/ScaleRov;
                        end
                        cfirstbest=c;
                        %  relative planner weights
                        for j=1:state
                            for k=1:agent
                                gammar(:,j,k)=gammagrid(:,agent+k);
                            end
                        end
                        % Compute the first best value functions without binding
                        % constraints
                        w=zeros(n,state,agent); dww=1;
                        wold=1/(1-beta)*utility(rho,cfirstbest(:,:,:));
                        while dww>0.00001
                            w(:,:,1)=utility(rho,cfirstbest(:,:,1))+beta*wold(:,:,1)*P';
                            w(:,:,2)=utility(rho,cfirstbest(:,:,2))+beta*wold(:,:,2)*P';
                            dww=max(max(max(abs(w-wold))));
                            wold=w;
                        end
                        wfirstbest=w;
                        % continuation utility as function of planner weight and state of the
                        % economy
                        ewfirstbest(:,:,1)=wfirstbest(:,:,1)*P';
                        ewfirstbest(:,:,2)=wfirstbest(:,:,2)*P';
                        
                        % ====================================================
                        % Start the loop with the enforcement constraints
                        % ====================================================n
                        ewsb=ewfirstbest;
                        wsb=wfirstbest;
                        options=optimset('Display','off');
                        eps=0.0000001;t=0;itgap=0; GAP=2; gap=2; gaphist=gap; gapincrease=0;
                        csolold=cfirstbest;
                        if beta<0.85
                            gapcrit=gapcritset;
                            maxitgap=maxitgapset;
                        else
                            gapcrit=gapcritset*1
                            maxitgap=maxitgapset;
                        end
                        
                        while GAP >=gapcrit && itgap<=maxitgap
                            rr=[];
                            rr1=[];
                            rr2=[];
                            rrnone=[];
                            rrboth=[];
                            nonconv=zeros(state,3)*nan;
                            itgap=itgap+1;
                            coalnonsust=nan(state,2);
                            count=0; 
                            for j=1:state
                                
                                %%%for each grid point, figure out which set of
                                %%%individuals from rest of village would want to come
                                %%%along with individual. 

                                
                                
                                index=[1 j];
                                % 1. constraint binding for agent 1
                                % (individual)
                                rr1=find(utility(rho,cfirstbest(:,j,1))+beta*ewsb(:,j,1)<waut(j,1) & utility(rho,cfirstbest(:,j,2))+beta*ewsb(:,j,2)>=waut(j,2)) ;
                                if vss>1 & coaldev==1 & size(rr1)>0
                                    rr1cd=[];
                                    for g=size(waut_check,4):-1:0 %%%check from the largest group down, 0 is an individual deviation
                                        
                                            if g>0
                                            villmatch_sub=NaN;
                                            
                                            for k=1:size(waut_check,3) 
                                            %check that
                                            %the subgroup under
                                            %consideration
                                            %is a stable group size
                                            n_type_sub=[];
                                            if n_type(1)>0 & n_type(2)==0
                                            n_type_sub(1)=sum(idio_dist_aut_all(j,1:3,k,g));
                                            n_type_sub(2)=0;
                                            end
                                            if n_type(1)>0 & n_type(2)>0
                                            n_type_sub(1)=sum(idio_dist_aut_all(j,1:3,k,g));
                                            n_type_sub(2)=sum(idio_dist_aut_all(j,4:6,k,g));
                                            end
                                            if n_type(1)==0 & n_type(2)>0
                                            n_type_sub(1)=0;
                                            n_type_sub(2)=sum(idio_dist_aut_all(j,1:3,k,g));
                                            end
                                            gg=0;
                                            while gg<=g
                                                gg=gg+1;
                                                rr=find(MATCHESKEEP(gg,:,g)==1);
                                                if length(rr)>0
                                                n_type_matches(1)=length(rr);
                                                else
                                                n_type_matches(1)=0;
                                                end
                                                rr=find(MATCHESKEEP(gg,:,g)==2);
                                                if length(rr)>0
                                                n_type_matches(2)=length(rr);
                                                else
                                                n_type_matches(2)=0;
                                                end
                                                if n_type_sub==n_type_matches
                                                    villmatch_sub=gg;
                                                    gg=g+5;
                                                end
                                                clear n_type_matches
                                            end
                                            if isnan(waut_check_ind(j,1,k,g))
                                            else
                                                if STABLEHETFINAL(kmatch,villmatch_sub,g)==1
                                                        coalition=[utility(rho,cfirstbest(rr1,j,1))+beta*ewsb(rr1,j,1)-waut_check_ind(j,1,k,g)];
                                                        for kk=1:size(waut_check,2)
                                                            if waut_check(j,kk,k,g)~=0
                                                                coalition=[coalition utility(rho,cfirstbest(rr1,j,2))+beta*ewsb(rr1,j,2)-waut_check(j,kk,k,g)];
                                                            end
                                                        end
                                                        %%%Step 1: find the grid points for which you are now trying to find a solution
                                                        
                                                        %%%%do the others want to come along?
                                                        rrcheck=coalition<0;
                                                        rrcheck=prod(rrcheck,2);
                                                        rrcheck=find(rrcheck>0);
                                                        
                                                        if size(rrcheck>0) 
                                                            
                                                            bind=1;
                                                            
                                                            coalitionsolve=[utility(rho,cfirstbest(:,j,1))+beta*ewsb(:,j,1)-waut_check_ind(j,1,k,g)];
                                                            for kk=1:size(waut_check,2)
                                                                if waut_check(j,kk,k,g)~=0
                                                                    coalitionsolve=[coalitionsolve utility(rho,cfirstbest(:,j,2))+beta*ewsb(:,j,2)-waut_check(j,kk,k,g)];
                                                                end
                                                            end
                                                            rrchecksolve=coalitionsolve<0;
                                                            rrchecksolve=prod(rrchecksolve,2);
                                                            rrchecksolve=find(rrchecksolve>0);
                                                            
                                                            rr=rrchecksolve(end)+1; 
                                                            if rr>length(gammagrid)
                                                                coalnonsust(j,1)=-2;
                                                                output=[nan,nan,-50];
                                                            else
                                                                
                                                                clear x0
                                                                
                                                                x0(1) = cfirstbest(rr,j,1);  %%redundant because we no longer need an initial guess
                                                                %T vectors of consumption and utilty that meet
                                                                %participation constraints
                                                                [x, output]=agent1maxrho_rov_06_Mar(x0,index,ewsb,[waut_check_ind(:,1,k,g) waut(:,2)],bind,Y,n,beta,rho,ScaleRov,cfirstbest,rr,coaldev);
                                                                
                                                                %T 29_Mar This keeps track of nonconvergence, but
                                                                %allows the programme to continue - sometimes it
                                                                %converges in a subsequent iteration
                                                            end
                                                            if output(3) <=0
                                                                nonconv(j,1)=output(3);
                                                                csol(rr1(rrcheck),j,1)=ones(size(rr1(rrcheck),1),1)*caut(j,1);
                                                                csol(rr1(rrcheck),j,2)=ones(size(rr1(rrcheck),1),1)*caut(j,2);
                                                                gammar(rr1(rrcheck),j,2)=(Y(j)-caut(j,1))/caut(j,1)/ScaleRov;
                                                                wsbnewsol(rr1(rrcheck),j,1)=ones(size(rr1(rrcheck),1),1)*waut(j,1);
                                                                wsbnewsol(rr1(rrcheck),j,2)=ones(size(rr1(rrcheck),1),1)*waut(j,2);
                                                                phisol(rr1(rrcheck),j,1)=gammagrid(rr1(rrcheck),4).^rho./gammar(rr1(rrcheck),j,2).^rho-1;
                                                                phisol(rr1(rrcheck),j,2:agent)=0	;
                                                            else
                                                                csol(rr1(rrcheck),j,1)=x(1);
                                                                gammar(rr1(rrcheck),j,2)=(Y(j)-x(1))/x(1)/ScaleRov;
                                                                csol(rr1(rrcheck),j,2)=(Y(j)-x(1))/ScaleRov;
                                                                wsbnewsol(rr1(rrcheck),j,1)=output(1);
                                                                wsbnewsol(rr1(rrcheck),j,2)=output(2);
                                                                phisol(rr1(rrcheck),j,1)=gammagrid(rr1(rrcheck),4).^rho./gammar(rr1(rrcheck),j,2).^rho-1;
                                                                phisol(rr1(rrcheck),j,2:agent)=0	;
                                                            end
                                                           
                                                        end
                                                end
                                                                  
                                                 
                                            end
                                            
                                            end
                                            rr1(rrcheck)=[];  %%%these grid points have been dealt with for coalitional deviation of size g
                                            rrnoneadd=rr1; %%%these are the grid points you still need to deal with for coalitional deviations.    
                                          
            
                                               
                                            elseif g==0 %%%only an individual deviation
                                                coalition=[utility(rho,cfirstbest(rr1,j,1))+beta*ewsb(rr1,j,1)-eeewaut(j,1)];
                                                
                                                rrcheck=coalition<0;
                                                rrcheck=prod(rrcheck,2);
                                                rrcheck=find(rrcheck>0);
                                                if size(rrcheck>0)
                                                    
                                                    bind=1;
                                                    coalitionsolve=[utility(rho,cfirstbest(:,j,1))+beta*ewsb(:,j,1)-eeewaut(j,1)];
                                                    %%%%do the others want to come along?
                                                    rrchecksolve=coalitionsolve<0;
                                                    rrchecksolve=prod(rrchecksolve,2);
                                                    rrchecksolve=find(rrchecksolve>0);
                                                    
                                                    rr=rrchecksolve(end)+1;
                                                    
                                                    if rr>length(gammagrid)
                                                        coalnonsust(j,1)=-2;
                                                        output=[nan,nan,-50];
                                                    else
                                                        clear x0
                                                        x0(1) = cfirstbest(rr,j,1);  %%redundant because we no longer need an initial guess
                                                        %T vectors of consumption and utilty that meet
                                                        %participation constraints
                                                        [x, output]=agent1maxrho_rov_06_Mar(x0,index,ewsb,[eeewaut],bind,Y,n,beta,rho,ScaleRov,cfirstbest,rr,coaldev);
                                                        
                                                        %T 29_Mar This keeps track of nonconvergence, but
                                                        %allows the programme to continue - sometimes it
                                                        %converges in a subsequent iteration
                                                    end
                                                    if output(3) <=0
                                                       nonconv(j,1)=output(3);
                                                            csol(rr1(rrcheck),j,1)=ones(size(rr1(rrcheck),1),1)*caut(j,1);
                                                        csol(rr1(rrcheck),j,2)=ones(size(rr1(rrcheck),1),1)*caut(j,2);
                                                        gammar(rr1(rrcheck),j,2)=(Y(j)-caut(j,1))/caut(j,1)/ScaleRov;
                                                        wsbnewsol(rr1(rrcheck),j,1)=ones(size(rr1(rrcheck),1),1)*waut(j,1);
                                                        wsbnewsol(rr1(rrcheck),j,2)=ones(size(rr1(rrcheck),1),1)*waut(j,2);
                                                        phisol(rr1(rrcheck),j,1)=gammagrid(rr1(rrcheck),4).^rho./gammar(rr1(rrcheck),j,2).^rho-1;
                                                        phisol(rr1(rrcheck),j,2:agent)=0	;
                                                    else
                                                        csol(rr1(rrcheck),j,1)=x(1);
                                                        gammar(rr1(rrcheck),j,2)=(Y(j)-x(1))/x(1)/ScaleRov;
                                                        csol(rr1(rrcheck),j,2)=(Y(j)-x(1))/ScaleRov;
                                                        wsbnewsol(rr1(rrcheck),j,1)=output(1);
                                                        wsbnewsol(rr1(rrcheck),j,2)=output(2);
                                                        phisol(rr1(rrcheck),j,1)=gammagrid(rr1(rrcheck),4).^rho./gammar(rr1(rrcheck),j,2).^rho-1;
                                                        phisol(rr1(rrcheck),j,2:agent)=0	;
                                                    end
                                                rr1(rrcheck)=[];  %%%these grid points have been dealt with
                                                rrnoneadd=rr1; %%these are additional gridpoints where noone has a binding constraint.
                                                end
                                                
                                            end
                                        
                                    end
                               
                                else
                                    %%%%in first group size iteration, i.e. for vss=1
                                    rr1cd=rr1;
                                    rr1id=[];
                                    rrnoneadd=[];
                                end

                                
                                if size(rr1cd)>0
                                    bind=1;
                                    rr=rr1cd(end)+1;
                                    if rr>length(gammagrid)
                                        coalnonsust(j,1)=-2;
                                        output=[nan,nan,-50];
                                        
                                    else
                                            clear x0
                                            x0(1) = cfirstbest(rr1cd(1),j,1);  %%redundant because we no longer need an initial guess
                                            %T vectors of consumption and utilty that meet
                                            %participation constraints
                                            [x, output]=agent1maxrho_rov_06_Mar(x0,index,ewsb,waut,bind,Y,n,beta,rho,ScaleRov,cfirstbest,rr,coaldev);
                                            
                                            %T 29_Mar This keeps track of nonconvergence, but
                                            %allows the programme to continue - sometimes it
                                            %converges in a subsequent iteration
                                    end
                                                              
                                    if output(3) <=0
                                        nonconv(j,1)=output(3);
                                        csol(rr1cd,j,1)=ones(size(rr1cd,1),1)*caut(j,1);
                                        csol(rr1cd,j,2)=ones(size(rr1cd,1),1)*caut(j,2);
                                        gammar(rr1cd,j,2)=(Y(j)-caut(j,1))/caut(j,1)/ScaleRov;
                                        wsbnewsol(rr1cd,j,1)=ones(size(rr1cd,1),1)*waut(j,1);
                                        wsbnewsol(rr1cd,j,2)=ones(size(rr1cd,1),1)*waut(j,2);
                                        phisol(rr1cd,j,1)=gammagrid(rr1cd,4).^rho./gammar(rr1cd,j,2).^rho-1;
                                        phisol(rr1cd,j,2:agent)=0	;
                                    else
                                        csol(rr1cd,j,1)=x(1);
                                        gammar(rr1cd,j,2)=(Y(j)-x(1))/x(1)/ScaleRov;
                                        csol(rr1cd,j,2)=(Y(j)-x(1))/ScaleRov;
                                        wsbnewsol(rr1cd,j,1)=output(1);
                                        wsbnewsol(rr1cd,j,2)=output(2);
                                        phisol(rr1cd,j,1)=gammagrid(rr1cd,4).^rho./gammar(rr1cd,j,2).^rho-1;
                                        phisol(rr1cd,j,2:agent)=0	;
                                    end
                                end
                                %%%%this becomes relevant for vss>1, these are all the left-overs that are not dealt with by joint deviations.
                                if size(rrnoneadd)>0
                                    
                                    csol(rrnoneadd,j,1:agent)=cfirstbest(rrnoneadd,j,1:agent);
                                    gammar(rrnoneadd,j,1:agent)=gammagrid(rrnoneadd,agent+1:2*agent);
                                    phisol(rrnoneadd,j,1:agent)=0;
                                    wsbnewsol(rrnoneadd,j,1)=utility(rho,cfirstbest(rrnoneadd,j,1))+beta*ewsb(rrnoneadd,j,1);
                                    wsbnewsol(rrnoneadd,j,2)=utility(rho,cfirstbest(rrnoneadd,j,2))+beta*ewsb(rrnoneadd,j,2);
                                end
                                
                                % 2. Constraint binding for agent 2 (RoV), 
                                rr2=find(utility(rho,cfirstbest(:,j,2))+beta*ewsb(:,j,2)<waut(j,2)) ;
                                
                                if size(rr2)>0
                                    bind=2; %%i.e. agent 2 has the binding constraint
                                    rr=find(utility(rho,cfirstbest(:,j,2))+beta*ewsb(:,j,2)>=waut(j,2));
                                    if isempty(rr)==1
                                        coalnonsust(j,2)=-2;
                                        output=[nan,nan,-50];
                                    else
                                        clear x0
                                        x0(1) = cfirstbest(rr(end),j,1); %%redundant because we no longer need an initial guess
                                        %%%%check whether 
                                        %%%%individual is happy with this allocation
                                        %%%%relative to autarky, given that rov gets
                                        %%%%waut(j,2) inside contract
                                        [x, output]=agent1maxrho_rov_06_Mar(x0,index,ewsb,[eeewaut(:,1) waut(:,2)],bind,Y,n,beta,rho,ScaleRov,cfirstbest,rr,coaldev);
                                        %%%%check that agent 1 does not have a binding
                                        %%%%constraint at this new solution, and note that
                                        %%%%the rov must never have less than waut(j,2),
                                        %%%%so if individual has a binding constraint that
                                        %%%%means rov would have to decrease its
                                        %%%%consumption, not possible
                                    end
                                    % output(3) equals 1 if constraints
                                    % are satisfied, but -2 if not
                                    if output(3)<=0
                                        nonconv(j,2)=output(3)*100;
                                        csol(rr2,j,1)=ones(size(size(rr2,1),1),1)*caut(j,1);
                                        gammar(rr2,j,2)=(Y(j)-caut(j,1))/caut(j,1)/ScaleRov;
                                        csol(rr2,j,2)=ones(size(size(rr2,1),1),1)*caut(j,2);
                                        wsbnewsol(rr2,j,1)=ones(size(rr2))*waut(j,1);
                                        wsbnewsol(rr2,j,2)=ones(size(rr2))*waut(j,2);
                                        phisol(rr2,j,2)=gammar(rr2,j,2).^rho./gammagrid(rr2,4).^rho-1;
                                        phisol(rr2,j,1)=0;
                                    else
                                        %%%%check that agent 1 (in conjunction with the rest of the village)
                                        %%%%truly does not have binding constraints
                                        %%%%at the proposed solution
                                     if vss>1 & coaldev==1 
                                       for g=size(waut_check,4):-1:1 %%%0 is an individual deviation and has already been dealt with above
                                        
                                            villmatch_sub=NaN;
                                            %check that
                                            %the subgroup you're looking at
                                            %is a stable group size
                                            for k=size(waut_check,3):-1:1   
                                            n_type_sub=[];
                                            if n_type(1)>0 & n_type(2)==0
                                            n_type_sub(1)=sum(idio_dist_aut_all(j,1:3,k,g));
                                            n_type_sub(2)=0;
                                            end
                                            if n_type(1)>0 & n_type(2)>0
                                            n_type_sub(1)=sum(idio_dist_aut_all(j,1:3,k,g));
                                            n_type_sub(2)=sum(idio_dist_aut_all(j,4:6,k,g));
                                            end
                                            if n_type(1)==0 & n_type(2)>0
                                            n_type_sub(1)=0;
                                            n_type_sub(2)=sum(idio_dist_aut_all(j,1:3,k,g));
                                            end
                                            gg=0;
                                            while gg<=g
                                                gg=gg+1;
                                                rr=find(MATCHESKEEP(gg,:,g)==1);
                                                if length(rr)>0
                                                n_type_matches(1)=length(rr);
                                                else
                                                n_type_matches(1)=0;
                                                end
                                                rr=find(MATCHESKEEP(gg,:,g)==2);
                                                if length(rr)>0
                                                n_type_matches(2)=length(rr);
                                                else
                                                n_type_matches(2)=0;
                                                end
                                                if n_type_sub==n_type_matches
                                                    villmatch_sub=gg;
                                                    gg=g+5;
                                                end
                                                clear n_type_matches
                                            end
                                            if isnan(waut_check_ind(j,1,k,g))
                                            else
                                                if STABLEHETFINAL(kmatch,villmatch_sub,g)==1
                                                        coalition=[output(1)-waut_check_ind(j,1,k,g)];
                                                        for kk=1:size(waut_check,2)
                                                            if waut_check(j,kk,k,g)~=0
                                                                coalition=[coalition output(2)-waut_check(j,kk,k,g)];
                                                            end
                                                        end
                                                        %%%Step 1: find the grid points for which you are now trying to find a solution
                                                        
                                                        %%%%do the others want to come along?
                                                        rrcheck=coalition<0;
                                                        rrcheck=prod(rrcheck,2);
                                                        rrcheck=find(rrcheck>0);
                                                        
                                                        if size(rrcheck>0) 
                                                         output(3)=-2;
                                                        end
                                                end
                                                end
                                           end 
  
                                       end
                                    end
                                        if output(3)>0
                                            csol(rr2,j,1)=x(1);
                                            gammar(rr2,j,2)=(Y(j)-x(1))/x(1)/ScaleRov;
                                            csol(rr2,j,2)=(Y(j)-x(1))/ScaleRov;
                                            wsbnewsol(rr2,j,1)=output(1);
                                            wsbnewsol(rr2,j,2)=output(2);
                                            phisol(rr2,j,2)=gammar(rr2,j,2).^rho./gammagrid(rr2,4).^rho-1;
                                            phisol(rr2,j,1)=0;
                                        else
                                            nonconv(j,2)=output(3)*100;
                                            csol(rr2,j,1)=ones(size(size(rr2,1),1),1)*caut(j,1);
                                            gammar(rr2,j,2)=(Y(j)-caut(j,1))/caut(j,1)/ScaleRov;
                                            csol(rr2,j,2)=ones(size(size(rr2,1),1),1)*caut(j,2);
                                            wsbnewsol(rr2,j,1)=ones(size(rr2))*waut(j,1);
                                            wsbnewsol(rr2,j,2)=ones(size(rr2))*waut(j,2);
                                            phisol(rr2,j,2)=gammar(rr2,j,2).^rho./gammagrid(rr2,4).^rho-1;
                                            phisol(rr2,j,1)=0;
                                        end
                                    end
                                end
                                
                                % 3. Constraint all slack
                                rrnone=find(utility(rho,cfirstbest(:,j,1))+beta*ewsb(:,j,1)>=waut(j,1) & utility(rho,cfirstbest(:,j,2))+beta*ewsb(:,j,2)>=waut(j,2)) ;
                                if size(rrnone)>0
                                    csol(rrnone,j,1:agent)=cfirstbest(rrnone,j,1:agent);
                                    gammar(rrnone,j,1:agent)=gammagrid(rrnone,agent+1:2*agent);
                                    phisol(rrnone,j,1:agent)=0;
                                   
                                    wsbnewsol(rrnone,j,1)=utility(rho,cfirstbest(rrnone,j,1))+beta*ewsb(rrnone,j,1);
                                    wsbnewsol(rrnone,j,2)=utility(rho,cfirstbest(rrnone,j,2))+beta*ewsb(rrnone,j,2);
                                end
                            end
                            gapold=max(norm(wsbnewsol(:,:,1)-wsb(:,:,1)),norm(wsbnewsol(:,:,2)-wsb(:,:,2)));
                            gap1=max(max(max(abs(invutility(beta,rho,wsbnewsol(:,:,1))-invutility(beta,rho,wsb(:,:,1)))./invutility(beta,rho,wsb(:,:,1)))),max(max(abs(invutility(beta,rho,wsbnewsol(:,:,2))-invutility(beta,rho,wsb(:,:,2)))./invutility(beta,rho,wsb(:,:,2)))))/(1-beta);
                            gap2=max(max(max(abs(csol(:,:,1)-csolold(:,:,1))./csol(:,:,1))),max(max(abs(csol(:,:,2)-csolold(:,:,2))./csol(:,:,2))));
                            gap=max(gap1,gap2);
                            gaphist(itgap+1)=gap;
                            GAP=max(gaphist(max(1,length(gaphist)-convlength):end))
                            disp('gap,gapold')
                            disp([gap,gapold])
                            nonconvind(Qbeta,Qrho,vss,incproc)=mean(nonconv(find(~isnan(nonconv))));
                            Coalnonsust(Qbeta,Qrho,vss,incproc)=mean(coalnonsust(find(~isnan(coalnonsust))));
                            
                         
                            if  (min(min(nonconv))<=0 ||min(min(coalnonsust))<0) & itgap>=3 & coaldev==1
                                itgap
                                disp('Autarky or coalition not sustainable')
                                [nonconv,S]
                                for j=1:state
                                    csol(:,j,1)=caut(j,1);
                                    csol(:,j,2)=caut(j,2);
                                end
                                break
                            end
                            csolold=csol;
                            wsb=wsbnewsol;
                            ewsb(:,:,1)=(P*wsb(:,:,1)')';
                            ewsb(:,:,2)=(P*wsb(:,:,2)')';
                        end
                        CRITGAP(vss)=gaphist(itgap-1);
                        if itgap>=maxitgap
                            nonconvind(Qbeta,Qrho,vss,incproc)=gap;
                        end
                        if gapincrease==1
                            nonconvind(Qbeta,Qrho,vss,incproc)=-gap;
                        end
                        
                        % =============================================================
                        % calculate values of outside option for coalitional deviations
                        % =============================================================
                        
                        if coaldev==1 && (nonconvind(Qbeta,Qrho,vss,incproc)<=0 || Coalnonsust(Qbeta,Qrho,vss,incproc)<=0) && itgap>=3 && vss>=1
                            % if vss is unsustainable, the outside option for vss+1 remains that
                            % for vss
                            waut_lag(:,1)=waut_lag(:,1);
                            waut_lag(:,2)=waut_lag(:,2);
                            P_lag=P_lag;
                            S_lag=S_lag;
                            % if the current insurance group is not sustainable, it
                            % cannot be the outside option for the next bigger group,
                            % so increase the 'lagind' indicator by one
                            lagind=lagind+1;
                            gap=-1;
                            STABLE(vss)=0;
                            STABLEHET(kmatch,villmatch,vss)=0;
                            STABLECHECK=[1 STABLE];
                            S_rov_KEEP(1:state/3,kmatch,villmatch,vss)=S_rov;
                            S_KEEP(1:state,:,kmatch,villmatch,vss)=S;
                            
                            
                            IDIO_DIST_ROV_KEEP(1:state/3,1:size(idio_dist,2),kmatch,villmatch,vss)=idio_dist;

                       
                        else 
                            lagind=1;
                            % when ins group of size vss is sustainable, its value
                            % becomes outside option for next bigger group size
                            % NOTE THAT EQUAL-SHARING IS NO LONGER
                            % NECESSARILY SITUATED AT n+1/2
                            dev=(gammagrid(:,4)-1).^2;
                            [CC,devindex]=min(dev);
                            for j=1:state
                                % outside option for next larger insurance group is
                                % value with equal weights vss ...
                                eewaut(j,1)=wsb(devindex,j,1);
                                eewaut(j,2)=wsb(devindex,j,2);
                                
                            end
                            waut_lag=eewaut;
                            S_lag=S;
                            P_lag=P;
                            S_rov_KEEP(1:state/3,kmatch,villmatch,vss)=S_rov;
                            S_KEEP(1:state,:,kmatch,villmatch,vss)=S;
                            C_KEEP(:,1:state,:,kmatch,villmatch,vss)=csol;
                            PHI_KEEP(:,1:state,:,kmatch,villmatch,vss)=phisol;
                            WAUT_LAG(1:state,1:2,kmatch,villmatch,vss)=waut_lag;
                            P_LAG(1:state,1:state,vss,kmatch,villmatch)=P_lag;
                            S_LAG(1:state,:,kmatch,villmatch,vss)=S_lag;
                            STABLE(vss)=1;
                            STABLEHET(kmatch,villmatch,vss)=1;
                            STABLECHECK=[1 STABLE]; %%%this includes stability for the individual as well, makes the loop above easier
                            IDIO_DIST_ROV_KEEP(1:state/3,1:size(idio_dist,2),kmatch,villmatch,vss)=idio_dist;
                        end
                        if itgap>=maxitgap
                            nonconvind(Qbeta,Qrho,vss,incproc)=gap;
                        end
                        
                 
                     
                
                
            end
            end
            %%%check consistency across village and individual when they
            %%%take different roles
            STABLEHETFINAL(1:2,1:size(MATCHES,1),vss)=STABLEHET(1:2,1:size(MATCHES,1),vss);
                for k=1:size(MATCHES,1)-1
                STABLEHETFINAL(1,k+1,vss)=STABLEHET(1,k+1,vss)*STABLEHET(2,k,vss);
                STABLEHETFINAL(2,k,vss)=STABLEHET(1,k+1,vss)*STABLEHET(2,k,vss);
                end
            end
            if Qbeta<=5
                 STABLEHETFINAL_BETA(1:2,1:5,1:4,Qbeta)=STABLEHETFINAL;
             else
                 STABLEHETFINAL_BETA(1:2,1:7,1:6,Qbeta)=STABLEHETFINAL;
             end
            
            
            %==========================================================
            % Simulate the economy
            % =========================================================
            % This lists all possible combinations of individual income
            % values (IDIO_DIST), and the aggregate states (pair of individual income and ROV-income) associated
            % with every individual income value
            
            if coaldev==1 
                %%%You find these by hand, but the rule is the largest
                %%%possible groups when the village is divided into
                %%%homogenous groups (i.e. low income and high income
                %%%groups) AND largest possible groups when the village is
                %%%divided into mixed groups (all the same mix) and no more than 1 difference in mix. 
            
                
            if beta==0.95
                sim_group_info(:,:,1)=[1 1; 2 2];
                sim_group_info(:,:,2)=[1 2; 2 1];
                sim_group_info(:,:,3)=[1 3; 2 2];
                vss_set_sim=[1 1 2];
            elseif beta==0.96
                sim_group_info(:,:,1)=[1 1; 2 2];
                sim_group_info(:,:,2)=[1 2; 2 1];
                sim_group_info(:,:,3)=[1 1; 2 3];
                sim_group_info(:,:,4)=[1 3; 2 2];
                vss_set_sim=[1 1 2 2];
            elseif beta==0.97
                sim_group_info(:,:,1)=[1 1; 2 2];
                sim_group_info(:,:,2)=[1 2; 2 1];
                sim_group_info(:,:,3)=[1 1; 2 3];
                sim_group_info(:,:,4)=[1 3; 2 2];
                sim_group_info(:,:,5)=[1 4; 2 3];
                vss_set_sim=[1 1 2 2 3];
            elseif beta==0.98
                sim_group_info(:,:,1)=[1 1; 2 2];
                sim_group_info(:,:,2)=[1 2; 2 1];
                sim_group_info(:,:,3)=[1 1; 2 3];
                sim_group_info(:,:,4)=[1 3; 2 2];        
                sim_group_info(:,:,5)=[1 1; 2 4];
                sim_group_info(:,:,6)=[1 4; 2 3];
                sim_group_info(:,:,7)=[1 5; 2 4];
                vss_set_sim=[1 1 2 2 3 3 4];
            
            elseif beta==0.9 | beta==0.91 | beta==0.92 | beta==0.93 | beta==0.94
                sim_group_info(:,:,1)=[1 1; 2 2];
                sim_group_info(:,:,2)=[1 2; 2 1];
                vss_set_sim=[1 1];
            end
            end
          for gg=1:length(vss_set_sim)
            vss=vss_set_sim(gg);  
            Nvill=ceil((vss_high_inddev+1)/(vss+1)); 
            cut=fix(Nvill/2);
            var_dyagglog=zeros(Nvill*(vss+1),vss_high);var_dyvillagglog=zeros(Nvill*(vss+1),vss_high);vardaggylog=zeros(Nvill*(vss+1),vss_high);betayc=zeros(2,Nvill*(vss+1));vardcshare=zeros(Nvill*(vss+1),vss_high);vardcshare_dylow=zeros(Nvill*(vss+1),vss_high);vardcshare_dyhigh=zeros(Nvill*(vss+1),vss_high);varclog=zeros(Nvill*(vss+1),vss_high);vardclog=zeros(Nvill*(vss+1),vss_high);varylog=zeros(Nvill*(vss+1),vss_high);vardylog=zeros(Nvill*(vss+1),vss_high);varclog_cond=zeros(Nvill*(vss+1),vss_high);varylog_cond=zeros(Nvill*(vss+1),vss_high);
                varc_yhigh_cond=zeros(Nvill*(vss+1),vss_high);varc_ylow_cond=zeros(Nvill*(vss+1),vss_high);vardclog_cond=zeros(Nvill*(vss+1),vss_high);vardylog_cond=zeros(Nvill*(vss+1),vss_high); varc_yhigh=zeros(Nvill*(vss+1),vss_high);
                vardc_dypos_cond=zeros(Nvill*(vss+1),vss_high);
                vardc_dyneg_cond=zeros(Nvill*(vss+1),vss_high); varc_yhigh=zeros(Nvill*(vss+1),vss_high); varc_ylow=zeros(Nvill*(vss+1),vss_high);vardc_dypos=zeros(Nvill*(vss+1),vss_high); vardc_dyneg=zeros(Nvill*(vss+1),vss_high);
                R2_dcdy=zeros(1,Nvill*(vss+1));betadcdy_agg=zeros(3,Nvill*(vss+1));   betadcdy=zeros(2,vss_high);  betadcdy_high=zeros(1,Nvill*(vss+1));betadcdy_low=zeros(1,Nvill*(vss+1));relcvar=zeros(Nvill*(vss+1),vss_high); reldcvar=zeros(Nvill*(vss+1),vss_high); reldc0var=zeros(Nvill*(vss+1),vss_high); relcvar_cond=zeros(Nvill*(vss+1),vss_high); reldcvar_cond=zeros(Nvill*(vss+1),vss_high);
                beta_cond=zeros(2,Nvill*(vss+1),vss_high); beta_ylow_cond=zeros(4,Nvill*(vss+1),vss_high); beta_yhigh_cond=zeros(4,Nvill*(vss+1),vss_high) ;
                beta_agg=zeros(2,Nvill*(vss+1),vss_high); beta_yhigh_agg=zeros(2,Nvill*(vss+1),vss_high); beta_ylow_agg=zeros(2,Nvill*(vss+1),vss_high);
                
                %=============================================
                % Start simulation
                % ============================================
                % number of insurance groups equals village size
                % vss_high+1 divided by size of ins group vss+1
                
                Npanel=Nvill*(vss+1);
                csim0=nan(Nvill,vss+1,T);
                income0=nan(Nvill,vss+1,T);
                income0ind=nan(Nvill,vss+1,T);
                phisim0=nan(Nvill,vss+1,T);
                Yvill=nan(1,T);
                Yagg=nan(Nvill,T);
                Ssim=nan(Nvill,vss+1,T);
                gammarfun=nan(Nvill,vss+1,T);
                phinew=nan(Nvill,vss+1,T);
                indSjj=nan(T,vss+1);
                disp('simulating') 
                disp('Reminder: now simulating for village,coaldev,incproc')
                disp([village coaldev incproc])
                disp('now simulating for parameters')
                disp([Qbeta,Qrho])
                disp('vss')
                disp(vss)
            for kk=1:2
                    kmatch=sim_group_info(kk,1,gg);
                    villmatch=sim_group_info(kk,2,gg);
                    
                    n_type(1)=length(find(MATCHESKEEP(villmatch,:,vss)==1)); %%%how many of each income type
                    n_type(2)=length(find(MATCHESKEEP(villmatch,:,vss)==2)); %%%%how many of each income type     
                    
                    n_type_all=n_type;
                    if kmatch==1
                        n_type_all(1)=n_type(1)+1;
                    elseif kmatch==2
                        n_type_all(2)=n_type(2)+1;
                    end
                    %%%do we need to split the village for simulation (yes,
                    %%%if we have homogeneous groups, no if we have
                    %%%heterogeneous groups
                    if n_type_all(1)>0 & n_type_all(2)>0 
                    if kk==1
                    composition_cut=[1:Nvill]
                    else
                    composition_cut=[];
                    end
                    else    
                    if kk==1
                    composition_cut=[1:cut];
                    elseif kk==2
                    composition_cut=[cut+1:Nvill];
                    end
                    end
                    
                    if n_type_all(1)>0 & n_type_all(2)>0
                       if kmatch==1
                       kmatch_aux=kmatch+1;
                       villmatch_aux=villmatch-1;
                       elseif kmatch==2
                       kmatch_aux=kmatch-1;
                       villmatch_aux=villmatch+1;
                       end
                       n_type_aux(1)=length(find(MATCHESKEEP(villmatch_aux,:,vss)==1)); %%%how many of each income type
                       n_type_aux(2)=length(find(MATCHESKEEP(villmatch_aux,:,vss)==2)); %%%%how many of each income type     
                   
                    end
                   
                    if rov_average==1
                        ScaleRov=vss;
                    else
                        ScaleRov=1;
                    end
                     
                    state=STATE(kmatch, villmatch,vss);
                    if n_type_all(1)>0 & n_type_all(2)>0
                    state_aux=STATE(kmatch_aux,villmatch_aux,vss);
                    end
                    
                    S=S_KEEP(1:state,:,kmatch,villmatch,vss);  
                    if n_type(1)>0 & n_type(2)>0
                        idio_dist=IDIO_DIST_ROV_KEEP(1:state/3,1:6,kmatch,villmatch,vss);
                    else
                        idio_dist=IDIO_DIST_ROV_KEEP(1:state/3,1:3,kmatch,villmatch,vss);
                    end
                    if n_type_all(1)>0 & n_type_all(2)>0
                    S_aux=S_KEEP(1:state_aux,:,kmatch_aux,villmatch_aux,vss); 
                    if n_type_aux(1)>0 & n_type_aux(2)>0
                    idio_dist_aux=IDIO_DIST_ROV_KEEP(1:state_aux/3,1:6,kmatch_aux,villmatch_aux,vss);
                    else
                    idio_dist_aux=IDIO_DIST_ROV_KEEP(1:state_aux/3,1:3,kmatch_aux,villmatch_aux,vss);    
                    end
                    end
                    
                    phisol=PHI_KEEP(:,1:state,:,kmatch,villmatch,vss);
                    csol=C_KEEP(:,1:state,:,kmatch,villmatch,vss);
                    if n_type_all(1)>0 & n_type_all(2)>0
                    phisol_aux=PHI_KEEP(:,1:state_aux,:,kmatch_aux,villmatch_aux,vss); 
                    csol_aux=C_KEEP(:,1:state_aux,:,kmatch_aux,villmatch_aux,vss);  
                    end
                    %%%make sure you get the right IDIO_DIST
                    if kmatch==1 & n_type(1)>0 & n_type(2)>0
                    IDIO_DIST=[ones(state/3,1)*[1 0 0 0 0 0] + idio_dist; ones(state/3,1)*[0 1 0 0 0 0] + idio_dist; ones(state/3,1)*[0 0 1 0 0 0]+idio_dist];
                    elseif kmatch==2 & n_type(1)>0 & n_type(2)>0
                    IDIO_DIST=[ones(state/3,1)*[0 0 0 1 0 0] + idio_dist; ones(state/3,1)*[0 0 0 0 1 0] + idio_dist; ones(state/3,1)*[0 0 0 0 0 1]+idio_dist];
                    elseif (kmatch==1 & n_type(1)>0 & n_type(2)==0) 
                    IDIO_DIST=[ones(state/3,1)*[1 0 0] + idio_dist; ones(state/3,1)*[0 1 0] + idio_dist; ones(state/3,1)*[0 0 1]+idio_dist];
                    elseif  (kmatch==2 & n_type(1)==0)
                    IDIO_DIST=[ones(state/3,1)*[1 0 0] + idio_dist; ones(state/3,1)*[0 1 0] + idio_dist; ones(state/3,1)*[0 0 1]+idio_dist];
                    elseif kmatch==1 & n_type(1)==0 & n_type(2)>0
                    idio_dist=[ones(state/3,1)*[0 0 0] idio_dist];
                    IDIO_DIST=[ones(state/3,1)*[1 0 0 0 0 0] + idio_dist; ones(state/3,1)*[0 1 0 0 0 0] + idio_dist; ones(state/3,1)*[0 0 1 0 0 0]+idio_dist];
                    elseif kmatch==2 & n_type(1)>0 & n_type(2)==0 
                    idio_dist=[idio_dist ones(state/3,1)*[0 0 0]];
                    IDIO_DIST=[ones(state/3,1)*[0 0 0 1 0 0] + idio_dist; ones(state/3,1)*[0 0 0 0 1 0] + idio_dist; ones(state/3,1)*[0 0 0 0 0 1]+idio_dist];
                    end
                    
                    if n_type_all(1)>0 & n_type_all(2)>0  %%%in this case, we need auxiliary variables
                        if kmatch_aux==1 & n_type_aux(1)>0 & n_type_aux(2)>0
                        IDIO_DIST_aux=[ones(state_aux/3,1)*[1 0 0 0 0 0] + idio_dist_aux; ones(state_aux/3,1)*[0 1 0 0 0 0] + idio_dist_aux; ones(state_aux/3,1)*[0 0 1 0 0 0]+idio_dist_aux];
                        elseif kmatch_aux==2 & n_type_aux(1)>0 & n_type_aux(2)>0
                        IDIO_DIST_aux=[ones(state_aux/3,1)*[0 0 0 1 0 0] + idio_dist_aux; ones(state_aux/3,1)*[0 0 0 0 1 0] + idio_dist_aux; ones(state_aux/3,1)*[0 0 0 0 0 1]+idio_dist_aux];
                        elseif (kmatch_aux==1 & n_type_aux(2)==0) | (kmatch_aux==2 & n_type_aux(1)==0) %%this can't happen for aux
                        IDIO_DIST=[ones(state_aux/3,1)*[1 0 0] + idio_dist_aux; ones(state_aux/3,1)*[0 1 0] + idio_dist_aux; ones(state_aux/3,1)*[0 0 1]+idio_dist_aux];
                        elseif kmatch_aux==1 & n_type_aux(1)==0 & n_type_aux(2)>0
                        idio_dist_aux=[ones(state_aux/3,1)*[0 0 0] idio_dist_aux];
                        IDIO_DIST_aux=[ones(state_aux/3,1)*[1 0 0 0 0 0] + idio_dist_aux; ones(state_aux/3,1)*[0 1 0 0 0 0] + idio_dist_aux; ones(state_aux/3,1)*[0 0 1 0 0 0]+idio_dist_aux];
                        elseif kmatch_aux==2 & n_type_aux(1)>0 & n_type_aux(2)==0
                        idio_dist_aux=[idio_dist_aux ones(state_aux/3,1)*[0 0 0]];
                        IDIO_DIST_aux=[ones(state_aux/3,1)*[0 0 0 1 0 0] + idio_dist_aux; ones(state_aux/3,1)*[0 0 0 0 1 0] + idio_dist_aux; ones(state_aux/3,1)*[0 0 0 0 0 1]+idio_dist_aux];
                        end
                    end
                        
                        
                        
                        
                        
                    distribution=size(IDIO_DIST,1);
               
                
                
                
                tic
                %initial state
                t=1;
                % initial values of income and consumption (set =
                % income)
                %stream = RandStream.getGlobalStream
                %reset(stream);
                %%%so now you always simulate twice%%%
                rng(3,'twister');
                shockinittemp=randi(3,Nvill*(vss+1),1);
                shockinit=reshape(shockinittemp,Nvill,vss+1);
                income0ind(:,:,1)=shockinit;
                
                
                incomevals_match1=incomevals_match(:,1);
                incomevals_match2=incomevals_match(:,2);
                income0(composition_cut,1:n_type_all(1),1)=incomevals_match1(income0ind(composition_cut,1:n_type_all(1),1));
                income0(composition_cut,n_type_all(1)+1:n_type_all(1)+n_type_all(2),1)=incomevals_match2(income0ind(composition_cut,n_type_all(1)+1:n_type_all(1)+n_type_all(2),1));
                csim0(composition_cut,:,1)=income0(composition_cut,:,1);
                 
              
                Yagg(composition_cut,t)=sum(income0(composition_cut,:,t),2); %this is aggregate income in the group, not in the village
                IDIO_DIST_sim=[];
                
                while t<T
                    t=t+1;
                    %%%%simulate next period income taking into account
                    %%%% persistence across time
                    rng(t,'twister');
                    shocktemp=rand(Nvill*(vss+1),1);
                    shock=reshape(shocktemp,Nvill,vss+1);
                    for k=1:n_type_all(1)
                        rr=(find(shock(composition_cut,k)<=cumQ_match(income0ind(composition_cut,k,t-1),1,1)));
                        income0ind(composition_cut(rr),k,t)=1;
                        rr=(find(cumQ_match(income0ind(composition_cut,k,t-1),1,1)<shock(composition_cut,k) & cumQ_match(income0ind(composition_cut,k,t-1),2,1)>=shock(composition_cut,k)));
                        income0ind(composition_cut(rr),k,t)=2;
                        rr=(find(cumQ_match(income0ind(composition_cut,k,t-1),2,1)<shock(composition_cut,k) & cumQ_match(income0ind(composition_cut,k,t-1),3,1)>=shock(composition_cut,k)));
                        income0ind(composition_cut(rr),k,t)=3;
                        income0(composition_cut,k,t)=incomevals_match(income0ind(composition_cut,k,t),1)  ;
                    end
                    for k=n_type_all(1)+1:n_type_all(1)+n_type_all(2)
                        rr=(find(shock(composition_cut,k)<=cumQ_match(income0ind(composition_cut,k,t-1),1,2)));
                        income0ind(composition_cut(rr),k,t)=1;
                        rr=(find(cumQ_match(income0ind(composition_cut,k,t-1),1,2)<shock(composition_cut,k) & cumQ_match(income0ind(composition_cut,k,t-1),2,2)>=shock(composition_cut,k)));
                        income0ind(composition_cut(rr),k,t)=2;
                        rr=(find(cumQ_match(income0ind(composition_cut,k,t-1),2,2)<shock(composition_cut,k) & cumQ_match(income0ind(composition_cut,k,t-1),3,2)>=shock(composition_cut,k)));
                        income0ind(composition_cut(rr),k,t)=3;
                        income0(composition_cut,k,t)=incomevals_match(income0ind(composition_cut,k,t),2)  ;
                    end
                    Yagg(composition_cut,t)=sum(income0(composition_cut,:,t),2);
                end
                
                
                    t=1;
                    while t<T
                        t=t+1;
                        gammarfun(composition_cut,:,t)=(Yagg(composition_cut,t-1)*ones(1,vss+1)-csim0(composition_cut,:,t-1))./csim0(composition_cut,:,t-1)/ScaleRov;
                          
                        for ivill=composition_cut
                            if (n_type_all(1)>0 & n_type_all(2)==0) | (n_type_all(1)==0 & n_type_all(2)>0) 
                            IDIO_DIST_sim=find(((IDIO_DIST(:,1,1)==(sum(income0ind(ivill,:,t)==1))).*(IDIO_DIST(:,2,1)==(sum(income0ind(ivill,:,t)==2))).*(IDIO_DIST(:,3,1)==(sum(income0ind(ivill,:,t)==3))))==1)';
                            else
                            IDIO_DIST_sim=find(((IDIO_DIST(:,1)==(sum(income0ind(ivill,1:n_type_all(1),t)==1))).*(IDIO_DIST(:,2,1)==(sum(income0ind(ivill,1:n_type_all(1),t)==2))).*(IDIO_DIST(:,3,1)==(sum(income0ind(ivill,1:n_type_all(1),t)==3))).*(IDIO_DIST(:,4,1)==(sum(income0ind(ivill,n_type_all(1)+1:n_type_all(1)+n_type_all(2),t)==1))).*(IDIO_DIST(:,5,1)==(sum(income0ind(ivill,n_type_all(1)+1:n_type_all(1)+n_type_all(2),t)==2))).*(IDIO_DIST(:,6,1)==(sum(income0ind(ivill,n_type_all(1)+1:n_type_all(1)+n_type_all(2),t)==3))))==1)';
                            IDIO_DIST_sim_aux=find(((IDIO_DIST_aux(:,1)==(sum(income0ind(ivill,1:n_type_all(1),t)==1))).*(IDIO_DIST_aux(:,2,1)==(sum(income0ind(ivill,1:n_type_all(1),t)==2))).*(IDIO_DIST_aux(:,3,1)==(sum(income0ind(ivill,1:n_type_all(1),t)==3))).*(IDIO_DIST_aux(:,4,1)==(sum(income0ind(ivill,n_type_all(1)+1:n_type_all(1)+n_type_all(2),t)==1))).*(IDIO_DIST_aux(:,5,1)==(sum(income0ind(ivill,n_type_all(1)+1:n_type_all(1)+n_type_all(2),t)==2))).*(IDIO_DIST_aux(:,6,1)==(sum(income0ind(ivill,n_type_all(1)+1:n_type_all(1)+n_type_all(2),t)==3))))==1)';
                            end
                            for k=1:vss+1
                                % choose for individual k the aggregate state
                                % that applies
                                if n_type_all(1)>0 & n_type_all(2)>0 
                                if k<=n_type_all(1)
                                indSjj(t,k)=IDIO_DIST_sim(S(IDIO_DIST_sim,1)==income0(ivill,k,t)); 
                                gammarfun(ivill,k,t)=min(gammarfun(ivill,k,t),GAMMAMAX(kmatch,villmatch,vss)); 
                                gammarfun(ivill,k,t)=max(gammarfun(ivill,k,t),GAMMAMIN(kmatch,villmatch,vss));
                      
                                % check that aggregate and individual income0s are
                                % consistent with the states caculated above
                                if abs(Yagg(ivill,t)-(S(indSjj(t,k),1)+ScaleRov*S(indSjj(t,k),2)))>0.000001
                                    stop
                                end
                                elseif k>n_type_all(1)    
                                indSjj(t,k)=IDIO_DIST_sim_aux(S_aux(IDIO_DIST_sim_aux,1)==income0(ivill,k,t));
                                gammarfun(ivill,k,t)=min(gammarfun(ivill,k,t),GAMMAMAX(kmatch_aux,villmatch_aux,vss)); 
                                gammarfun(ivill,k,t)=max(gammarfun(ivill,k,t),GAMMAMIN(kmatch_aux,villmatch_aux,vss));
                      
                                if abs(Yagg(ivill,t)-(S_aux(indSjj(t,k),1)+ScaleRov*S_aux(indSjj(t,k),2)))>0.000001
                                    stop
                                end
                                end
                                
                                else
                                indSjj(t,k)=IDIO_DIST_sim(S(IDIO_DIST_sim,1)==income0(ivill,k,t));
                                gammarfun(ivill,k,t)=min(gammarfun(ivill,k,t),GAMMAMAX(kmatch,villmatch,vss)); 
                                gammarfun(ivill,k,t)=max(gammarfun(ivill,k,t),GAMMAMIN(kmatch,villmatch,vss));
                      
                                if abs(Yagg(ivill,t)-(S(indSjj(t,k),1)+ScaleRov*S(indSjj(t,k),2)))>0.000001
                                    stop
                                end
                                end
                                
                            end
                            Ssim(ivill,:,t)=indSjj(t,:);
                        end
                        for ivill=composition_cut
                        for k=1:vss+1
                           if n_type_all(1)>0 & n_type_all(2)>0
                           if k<=n_type_all(1)
                           gammagrid=[linspace(min(S(:,1)),1.2*max(S(:,1)),m)',linspace(1.2*max(S(:,2)),min(S(:,2)),m)',ones(m,1),linspace(1.2*max(S(:,2)),min(S(:,2)),m)'./linspace(min(S(:,1)),1.2*max(S(:,1)),m)'];
                           phinew(ivill,k,t)=interp1(gammagrid(:,4),phisol(:,Ssim(ivill,k,t),1),gammarfun(ivill,k,t));
                     
                           elseif k> n_type_all(1)
                           gammagrid=[linspace(min(S_aux(:,1)),1.2*max(S_aux(:,1)),m)',linspace(1.2*max(S_aux(:,2)),min(S_aux(:,2)),m)',ones(m,1),linspace(1.2*max(S_aux(:,2)),min(S_aux(:,2)),m)'./linspace(min(S_aux(:,1)),1.2*max(S_aux(:,1)),m)'];
                        
                           phinew(ivill,k,t)=interp1(gammagrid(:,4),phisol_aux(:,Ssim(ivill,k,t),1),gammarfun(ivill,k,t));
                         
                           end
                           else
                           gammagrid=[linspace(min(S(:,1)),1.2*max(S(:,1)),m)',linspace(1.2*max(S(:,2)),min(S(:,2)),m)',ones(m,1),linspace(1.2*max(S(:,2)),min(S(:,2)),m)'./linspace(min(S(:,1)),1.2*max(S(:,1)),m)'];
                            
                           phinew(ivill,k,t)=interp1(gammagrid(:,4),phisol(:,Ssim(ivill,k,t),1),gammarfun(ivill,k,t));
                          
                           end
                        end
                        end
                        gammaraux=((1+phinew(composition_cut,:,t))./(1+phinew(composition_cut,1,t)*ones(1,vss+1))).^(1/rho).*csim0(composition_cut,:,t-1)./(csim0(composition_cut,1,t-1)*ones(1,vss+1));

                        csim0(composition_cut,:,t)=(gammaraux./(sum(gammaraux,2)*ones(1,vss+1))).*(Yagg(composition_cut,t)*ones(1,vss+1));
                       
                    end
            end
                    
                    
                    AA=[];
                    BB=[];
                    CC=[];
                    for t=1:T
                        AA=[AA reshape(csim0(:,:,t)',Nvill*(vss+1),1)];
                        BB=[BB reshape(income0(:,:,t)',Nvill*(vss+1),1)];
                        CC=[CC reshape(income0ind(:,:,t)',Nvill*(vss+1),1)];
                    end
                    csim0keep=csim0;
                    csim0=AA;
                    income0ind=CC;
                    income0=BB;
                    
               
                    
               
                disp('simulating takes')
                tsim(Qbeta,Qrho)=toc;            
                disp(tsim(Qbeta,Qrho))
                
                % =========================================================
                % calculate moments of consumption and income growth
                % =========================================================
                csim =[]; cshare_raw =[]; dcshare_raw =[]; cc =[]; ccc =[]; yy =[]; aggyy =[]; aggvillincome=[]; yyy_cond  =[]; dc_log_cond =[]; dy_log_cond =[]; ccc_cond=[];
                rng(3, 'twister'); cerr=randn(size(csim0)); logcsim0_groupmeandev=[]; logcsim_groupmeandev=[]; logincome_groupmeandev=[];
                rng(4, 'twister'); incerr=randn(size(csim0));   Momentsest1=nan(35,N_cerr); ind=[];
                % loop over size of cons measurement error
                for kkk=1:size(csim0,1)
                    tempvar(kkk)=var(log(csim0(kkk,:)));
                end
                Tempvar=mean(tempvar);
                for ic=1:N_cerr+1
                    if N_cerr>0
                        csim=exp(log(csim0)+((ic-1)/N_cerr*maxcerr/(1-(ic-1)/N_cerr*maxcerr)*Tempvar)^0.5*cerr);
                    else
                        csim=exp(log(csim0));
                    end
                    logcsim=log(csim);
                    logcsim0=log(csim0);
                    logcsim0_meandev=logcsim0-ones(size(csim0,1),1)*mean(logcsim0,1);
                    logcsim_meandev=logcsim-ones(size(csim0,1),1)*mean(logcsim,1);
		    
                    % add inc meas error back in
                    income=exp(log(income0)+varnu(incproc)^0.5*incerr);
                    logincome=log(income);
                    logincome0=log(income0);
                    logincome_meandev=logincome-ones(size(csim0,1),1)*mean(logincome,1);
                    income0_meandev=income0-ones(size(csim0,1),1)*mean(income0,1);
                    for ivill=1:Nvill
                        aggvillincome((ivill-1)*(vss+1)+1:ivill*(vss+1),1:T)=ones(vss+1,1)*sum(income((ivill-1)*(vss+1)+1:ivill*(vss+1),:),1);
                        logaggvillincome=log(aggvillincome);
                        logcsim0_groupmeandev((ivill-1)*(vss+1)+1:ivill*(vss+1),1:T)=logcsim0((ivill-1)*(vss+1)+1:ivill*(vss+1),:)-ones(vss+1,1)*mean(logcsim0((ivill-1)*(vss+1)+1:ivill*(vss+1),:),1);
                        logcsim_groupmeandev((ivill-1)*(vss+1)+1:ivill*(vss+1),1:T)=logcsim((ivill-1)*(vss+1)+1:ivill*(vss+1),:)-ones(vss+1,1)*mean(logcsim((ivill-1)*(vss+1)+1:ivill*(vss+1),:),1);
                        logincome_groupmeandev((ivill-1)*(vss+1)+1:ivill*(vss+1),1:T)=logincome((ivill-1)*(vss+1)+1:ivill*(vss+1),1:T)-ones(vss+1,1)*mean(logincome((ivill-1)*(vss+1)+1:ivill*(vss+1),:),1);
                    end
                    aggincome=sum(income,1);
                    % this chooses log or level consumption for
                    % moments
                    if momentclevinlog==1
                        ccc=logcsim(:,Tinit:T);
                        ccc_cond=logcsim_meandev(:,Tinit:T);
                        ccc_groupcond=logcsim_groupmeandev(:,Tinit:T);
                    else
                        ccc=(csim(:,Tinit:T));
                        ccc_cond=(csim_meandev(:,Tinit:T)); %NB this doesn't actually exist...
                    end
                    yyy=logincome(:,Tinit:T);
                    yyy_cond=logincome_meandev(:,Tinit:T);
                    yyy_groupcond=logincome_groupmeandev(:,Tinit:T);
                    cc=logcsim(:,Tinit:T)-logcsim(:,Tinit-1:T-1);
	             var_dc_t_fe(vss)=var(mean(cc,1));
                    cc_cond=logcsim_meandev(:,Tinit:T)-logcsim_meandev(:,Tinit-1:T-1);
                    cc_groupcond=logcsim_groupmeandev(:,Tinit:T)-logcsim_groupmeandev(:,Tinit-1:T-1);
                    yy=logincome(:,Tinit:T)-logincome(:,Tinit-1:T-1);
                    yy_cond=logincome_meandev(:,Tinit:T)-logincome_meandev(:,Tinit-1:T-1);
                    yy_groupcond=logincome_meandev(:,Tinit:T)-logincome_meandev(:,Tinit-1:T-1);
                    aggyy=log(sum(income(:,Tinit:T),1))-log(sum(income(:,Tinit-1:T-1),1));
                    aggvillyy=log(aggvillincome(:,Tinit:T))-log(aggvillincome(:,Tinit-1:T-1));
                    vardaggylog(1,vss)=var(log(aggincome(1,Tinit:T))-log(aggincome(1,Tinit-1:T-1)));
                    
                    % calculate moments for each individiual
                    for ind=1:size(csim,1)
                        dc_log_cond(ind,:)=(ccc_cond(ind,2:end))-(ccc_cond(ind,1:end-1));
                        dy_log_cond(ind,:)=(yyy_cond(ind,2:end))-(yyy_cond(ind,1:end-1));
                        dc_log_groupcond(ind,:)=(ccc_groupcond(ind,2:end))-(ccc_groupcond(ind,1:end-1));
                        dy_log_groupcond(ind,:)=(yyy_groupcond(ind,2:end))-(yyy_groupcond(ind,1:end-1));
                        
                        % conditioning on income increases, and income
                        % decreases
                        % NB: we allocate no change to decrease
                        rrhc =[]; rrlc =[]; rrl =[]; rrh=[];
                        rrhc=find(yy(ind,:)>0);
                        rrlc=find(yy(ind,:)<=0);
                        rrhc_cond=find(yy_cond(ind,:)>0);
                        rrlc_cond=find(yy_cond(ind,:)<=0);
                        rrhc_groupcond=find(yy_groupcond(ind,:)>0);
                        rrlc_groupcond=find(yy_groupcond(ind,:)<=0);
                        vardclog(ind,vss)=var(cc(ind,:));
                        vardylog(ind,vss)=var(yy(ind,:));
                        vardc_dypos(ind,vss)=var(cc(ind,rrhc));
                        vardc_dyneg(ind,vss)=var(cc(ind,rrlc));
                        
                        %%%conditional on village income
                        vardclog_cond(ind,vss)=var(dc_log_cond(ind,:));
                        vardylog_cond(ind,vss)=var(dy_log_cond(ind,:));
                        vardc_dypos_cond(ind,vss)=var(dc_log_cond(ind,dy_log_cond(ind,:)>0));
                        vardc_dyneg_cond(ind,vss)=var(dc_log_cond(ind,dy_log_cond(ind,:)<=0));
                        
                        %%%conditional on group income
                        vardclog_groupcond(ind,vss)=var(dc_log_groupcond(ind,:));
                        vardylog_groupcond(ind,vss)=var(dy_log_groupcond(ind,:));
                        vardc_dypos_groupcond(ind,vss)=var(dc_log_groupcond(ind,dy_log_groupcond(ind,:)>0));
                        vardc_dyneg_groupcond(ind,vss)=var(dc_log_groupcond(ind,dy_log_groupcond(ind,:)<=0));
                      
                        vardc_dypos(ind,vss)=var(cc(ind,rrhc));
                        vardc_dyneg(ind,vss)=var(cc(ind,rrlc));
                        
                        %%% Relative variances
                        reldcvar(ind,vss)=var(cc(ind,rrhc))-var(cc(ind,rrlc));
                        reldcvar_cond(ind,vss)=vardc_dypos_cond(ind,vss)-vardc_dyneg_cond(ind,vss);
                        reldcvar_groupcond(ind,vss)=vardc_dypos_groupcond(ind,vss)-vardc_dyneg_groupcond(ind,vss);
                        
                        %%% THE BETA COEFFICIENTS
                        %1) Unconditional
                        %risk-sharing whole sample
                        covtemp=cov(cc(ind,:),aggyy(1,:));
                        betadcdy_agg(ind,vss)=covtemp(2,1)/var(aggyy(1,:));
                        covtemp=cov(cc(ind,:),yy(ind,:));
                        betadcdy(ind,vss)=covtemp(2,1)/(var(yy(ind,:)));
                        [betadcdytemp,temp1,rtemp,temp2,statstemp]=regress(cc(ind,:)',[ones(1,size(yy,2))',aggyy(1,:)',yy(ind,:)']);
                        rrstat(ind,vss)=rtemp(2);
                        betadcdyalt(ind,vss)=betadcdytemp(2);
                        RR2_dcdy(ind,vss)=statstemp(1);
                        clear rtemp betadcdytemp statstemp
			[betadcdytemp,temp1,rtemp,temp2,statstemp]=regress(cc_cond(ind,:)',[ones(1,size(yy,2))',yy_cond(ind,:)']);
                        rrstat(ind,vss)=rtemp(2);
                        betadcdyalt1(ind,vss)=betadcdytemp(2);
                        RR2_dcdy_cond(ind,vss)=statstemp(1);
                        clear rtemp betadcdytemp statstemp

                        % for those with income gains versus
                        %those with income losses
                        YY=cc(ind,:)';
                        yy_hc(ind,:)=zeros(1,size(cc,2));
                        yy_hc(ind,rrhc)=yy(ind,rrhc);
                        oneone=ones(1,size(cc,2));
                        ones_hc(ind,:)=zeros(1,size(cc,2));
                        ones_hc(ind,rrhc)=oneone(rrhc);
                        XX=[ones(1,size(cc,2))',ones_hc(ind,:)',yy(ind,:)',yy_hc(ind,:)'];
                        beta_yhigh(:,ind,vss)=regress(YY,XX);
                        betadcdy_low(ind,vss)=beta_yhigh(3,ind,vss);
                        betadcdy_high(ind,vss)=beta_yhigh(4,ind,vss)+beta_yhigh(3,ind,vss);
                        
                        %2) conditional on village income
                        covtemp=[];
                        X=[ones(size(yy_cond(ind,:),2),1),yy_cond(ind,:)'];%
                        YY=cc_cond(ind,:)';
                        beta_cond(:,ind,vss)=inv(X'*X)*(X'*YY);
                        X=[ones(size(aggyy(1,:),2),1),aggyy'];%
                        YY=cc(ind,:)';
                        beta_agg(:,ind,vss)=inv(X'*X)*(X'*YY);
                        %for those with income gains versus
                        %those with income losses
                        YY_cond=cc_cond(ind,:)';
                        yy_hc_cond(ind,:)=zeros(1,size(cc_cond,2));
                        yy_hc_cond(ind,rrhc_cond)=yy_cond(ind,rrhc_cond);
                        oneone=ones(1,size(cc_cond,2));
                        ones_hc_cond(ind,:)=zeros(1,size(cc_cond,2));
                        ones_hc_cond(ind,rrhc_cond)=oneone(rrhc_cond);
                        XX_cond=[ones(1,size(cc_cond,2))',ones_hc_cond(ind,:)',yy_cond(ind,:)',yy_hc_cond(ind,:)'];
                        beta_yhigh_cond(:,ind,vss)=regress(YY_cond,XX_cond);

                        %3) conditional on group income
                        covtemp=[];
                        X=[ones(size(yy_groupcond(ind,:),2),1),yy_groupcond(ind,:)'];%
                        YY=cc_groupcond(ind,:)';
                        beta_groupcond(:,ind,vss)=inv(X'*X)*(X'*YY);
                        X=[ones(size(aggvillyy(ind,:),2),1),aggvillyy(ind,:)'];%
                        YY=cc(ind,:)';
                        beta_groupagg(:,ind,vss)=inv(X'*X)*(X'*YY);
                        %risk-sharing for those with income gains versus
                        %those with income losses
                        YY_groupcond=cc_groupcond(ind,:)';
                        yy_hc_groupcond(ind,:)=zeros(1,size(cc_groupcond,2));
                        yy_hc_groupcond(ind,rrhc_groupcond)=yy_groupcond(ind,rrhc_groupcond);
                        oneone=ones(1,size(cc_groupcond,2));
                        ones_hc_groupcond(ind,:)=zeros(1,size(cc_groupcond,2));
                        ones_hc_groupcond(ind,rrhc_groupcond)=oneone(rrhc_groupcond);
                        XX_groupcond=[ones(1,size(cc_groupcond,2))',ones_hc_groupcond(ind,:)',yy_groupcond(ind,:)',yy_hc_groupcond(ind,:)'];
                        beta_yhigh_groupcond(:,ind,vss)=regress(YY_groupcond,XX_groupcond);
                    end
                    var_dyagglog(1,vss)=var(aggyy);
                    R2_dcdy=mean(RR2_dcdy_cond(1:size(csim,1),vss));
                    
                    % =====================
                    % Matrix of moments
                    % =====================                    
                    % col 1-5 are parameters
                    Momentsest1(1:5,ic)=[vss,beta,rho,rho1vec(incproc),varnu(incproc)]';
                    % 6-11 relative variance dc/dy, dc(dy>0)/dc(dy<0);  unconditional (6:7), conditional on village income
                    % (8:9), conditonal on group income (10:11)
                    Momentsest1(6:11,ic)=[mean(vardclog(1:size(csim,1),vss))/mean(vardylog(1:size(csim,1),vss)), mean(reldcvar(1:size(csim,1),vss))...
                    mean(vardclog_cond(1:size(csim,1),vss))/mean(vardylog_cond(1:size(csim,1),vss)), mean(reldcvar_cond(1:size(csim,1),vss))...
                    mean(vardclog_groupcond(1:size(csim,1),vss))/mean(vardylog_groupcond(1:size(csim,1),vss)), mean(reldcvar_groupcond(1:size(csim,1),vss))]';
                    % 12-14: beta, beta high, beta low unconditional
                    Momentsest1(12:14,ic)=[mean(betadcdy(1:size(csim,1),vss)),mean(betadcdy_high(1:size(csim,1),vss)),mean(betadcdy_low(1:size(csim,1),vss))]';
                    % 13_16: beta, beta high, beta low conditional on
                    % village income
                    Momentsest1(15:17,ic)=[mean(beta_cond(2,1:size(csim,1),vss),2),mean(beta_yhigh_cond(3,1:size(csim,1),vss),2)+mean(beta_yhigh_cond(4,1:size(csim,1),vss),2),mean(beta_yhigh_cond(3,1:size(csim,1),vss),2)]';
                    % 18:20: beta, beta high, beta low conditional on
                    % group income
                    Momentsest1(18:20,ic)=[mean(beta_groupcond(2,1:size(csim,1),vss),2),mean(beta_yhigh_groupcond(3,1:size(csim,1),vss),2)+mean(beta_yhigh_groupcond(4,1:size(csim,1),vss),2),mean(beta_yhigh_groupcond(3,1:size(csim,1),vss),2)]';
                    %21:22 are variances of dclog, dylog
                    Momentsest1(21:22,ic)=[mean(vardclog(1:size(csim,1),vss)),mean(vardylog(1:size(csim,1),vss))];
                    %23:28, relative variance dc/dy (dy>0) and var dc/dy(dy<0);
                    %unconditinoal 23:24, conditional on village income 25:26,
                    %conditional on group income 27:28
                    Momentsest1(23:24,ic)=[mean(vardc_dypos(1:size(csim,1),vss));mean( vardc_dyneg(1:size(csim,1),vss))];
                    Momentsest1(25:26,ic)=[mean(vardc_dypos_cond(1:size(csim,1),vss));mean(  vardc_dyneg_cond(1:size(csim,1),vss))];
                    Momentsest1(27:28,ic)=[mean(vardc_dypos_groupcond(1:size(csim,1),vss));mean( vardc_dyneg_groupcond(1:size(csim,1),vss))];
                    Momentsest1(30,ic)=[mean( vardylog_cond(1:size(csim,1),vss))];
                    Momentsest1(31,ic)=mean(RR2_dcdy(1:size(csim,1),vss));
                    Momentsest1(32,ic)=mean(vardclog(1:size(csim,1),vss));
                    Momentsest1(33,ic)=mean(vardclog_cond(1:size(csim,1),vss));
 		   Momentsest1(34,ic)=var_dc_t_fe(vss);
 		    Momentsest1(35,ic)=mean(RR2_dcdy_cond(1:size(csim,1),vss));
                    %%% and some other model-dependent measures
                    if coaldev<2
                        Momentsest1(29,ic)=CRITGAP(vss);
                    elseif coaldev==2
                        Momentsest1(29,ic)=phi;
                    elseif coaldev==3
                        Momentsest1(29,ic)=nan;
                    end
                end
                Momentsest(:,:,Qbeta,Qrho,vss,gg)=[Momentsest1];
                toc
           
        end
              
            

                cd(outputpath)
                save(filetosaveest,'Momentsest');
                cd(programpath)
        end
    end
  
    
    
%%%%%%%%%%%%%%%%%%%%%%%%%%%Momentsexplore 
 analytic=1
 outlier=0
  relvarincscale=1; 
 if Tessa1Tobi2==1
    outputpath='C:\Users\tbold\Dropbox\Tessa and Tobi\Documents\Latex';
    modeldatapath='C:\Users\tbold\Dropbox\Tessa and Tobi\JEEA Replication\Programmes\Main Results and Robustness';
    datapath='C:\Users\tbold\Dropbox\Tessa and Tobi\JEEA Replication/Data/Output';
    datapathvil1='C:\Users\tbold\Dropbox\Tessa and Tobi\JEEA Replication\Data/Output';
    addpath('C:\Users\tbold\Dropbox\Tessa and Tobi\JEEA Replication\Programmes');
 else
    outputpath='C:\Users\tbroer\Dropbox\Research\Current_Projects\tessa\Documents\Latex';
    modeldatapath='C:\Users\tbroer\Dropbox\Research\Current_Projects\tessa\JEEA Replication\Programmes\Main Results and Robustness';
    datapath='C:\Users\tbroer\Dropbox\Research\Current_Projects\tessa\JEEA Replication/Data/Output';
    datapathvil1='C:\Users\tbroer\Dropbox\Research\Current_Projects\tessa\JEEA Replication\Data/Output';
    addpath('C:\Users\tbroer\Dropbox\Research\Current_Projects\tessa\JEEA Replication\Programmes');
 end   
    datafile=['icrisatmoments' int2str(outlier)];
    if village==1
        cd(datapathvil1)
    else
        cd(datapath)
    end
    eval(['load ' datafile '.out'])
    eval(['datamoments=' datafile ';']);
    eval(['clear ' datafile])
    
    
    Nom=4;
    
    bootstrapvar=datamoments(1:Nom,4+(village-1)*Nom+1:4+(village)*Nom);
    analyticvar=datamoments(1:Nom,16+(village-1)*Nom+1:16+(village)*Nom);
    analyticmoments=datamoments(1:Nom,1:3);
    
    analyticrawmoments=datamoments(1:Nom,28+1:28+3);
    bootstraprawvar=datamoments(1:Nom,28+4+(village-1)*Nom+1:28+4+(village)*Nom);
    analyticrawvar=datamoments(1:Nom,28+16+(village-1)*Nom+1:28+16+(village)*Nom);
    %%%here you replace the exact standard errors (the variances are
    %%%bootstrapped in STATA, the coefficients are delta-method
    
    VCV2data=bootstrapvar;
    VCV2rawdata=bootstraprawvar;
    
    if analytic==1
        VCV2data(2,2)=analyticvar(2,2);
        VCV2data(4,4)=analyticvar(4,4);
        VCV2rawdata(2,2)=analyticrawvar(2,2);
        VCV2rawdata(4,4)=analyticrawvar(4,4);
    end
        
   
Momvec2data(1,1:4)=[0.30 0.21 -0.08 -0.41];
Vardylog_cond=[0.336;0.162;0.268];
        Momentssum=[];
      
            %By hand: for homogeneous groups, for Qbeta=1:6, vss=1, LL,HH
            %for Qbeta=7:8, vss=2, LLL HHH, Qbeta=9, vss=3, LLLL HHHH
            %are the largest homogeneous groups, these are selected from
            %the simulations by the index below
            index_homogeneous=[1 1 1; 2 1 1; 3 1 1; 4 1 1; 5 1 1; 6 1 1; 7 2 3; 8 2 3; 9 3 5];
            index_heterogeneous=[1 1 2; 2 1 2; 3 1 2; 4 1 2; 5 1 2; 6 1 2; 6 2 3; 7 1 2; 7 2 4; 8 1 2; 8 2 4; 8 3 5; 9 1 2; 9 2 4; 9 3 6; 9 4 7];  
            distribution_heterogeneous=[ 1 1; 1 1; 1 1; 1 1; 1 1; 1 1; 1 2; 1 1; 1 2; 1 1; 1 2; 1 3; 1 1; 1 2; 1 3; 1 4];
            homogeneous=0;  % set to 0 for estimating homogeneous groups and 1 for mixed groups     
          if homogeneous==1
              index=index_homogeneous;
          elseif homogeneous==0
              index=index_heterogeneous;
          end
          clear Momentstable W4 W4d
    
          for kk=1:size(index,1)
          inddelta=index(kk,1);
          indrho=1;
          vss=index(kk,2);
          incproc=1;
          j=1
          gg=index(kk,3);
          Momentstable(1:3,kk)=Momentsest(1:3,j,inddelta,indrho,vss,gg);
                
          Momentstable(4:7,kk)=[Momentsest([8 15],j,inddelta,indrho,vss,gg) ; Momentsest(9,j,inddelta,indrho,vss,gg)/Vardylog_cond(village) ;(Momentsest([16],j,inddelta,indrho,vss,gg)-Momentsest([17],j,inddelta,indrho,vss,gg))];
          Momentstable(8,kk)=(Momentsest(21,j,inddelta,indrho,vss,gg)-Momentsest(21,1,inddelta,indrho,vss,gg))./Momentsest(21,j,inddelta,indrho,vss,gg);%2*maxcerr*(j-1)/(size(moments,2)-1)/moments(28,j);
          Momentstable(9,kk)=2*Momentsest(5,j,inddelta,indrho,vss,gg)./vardy(village);
          Momentstable(1,kk)=Momentstable(1,kk)+1;
          MOMENTSDIAGNOSTICS(1:9,kk)=Momentstable(:,kk);
                
          MOMENTSDIAGNOSTICS(10,kk)=Momentsest(10,kk);
                
         momentsdev=Momentstable(4:7,kk)- Momvec2data(1,1:4)';
                            %%%2,4, d:diagonal
                            
                            
                            
          W4(kk)=momentsdev'*inv(VCV2data(1:4,1:4))*momentsdev;
          W4d(kk)=momentsdev'*inv(diag(diag(VCV2data(1:4,1:4))))*momentsdev;
                            
                               
                            
         end
                       
      
                
    [Momentstable; W4d]
    if homogeneous==1
    %these are the results bearing in mind that there must not be a difference greater 1 for income type    
    [Momentstable(:,9); W4d(9)]    
    else
    [Momentstable(:,9); W4d(9)]      
    end